LANGLEY RESEARCH CENTER 


1176 00501 6226 



NASA Contractor Report 16B816 




NASA-CR-165816 

19840021082 


A NUMERICAL SIMULATION OF THE DISPERSAL 
OF AERIAL SPRAYS 


M. B. Bragg 


THE OHIO STATE UNIVERSITY 

Aeronautical and Astronautical Research Laboratory 
Columbus, Ohio 43220 



NHi 


L-. , --- - H r-.^TEn: 

NASA Purchase Order L-7965B , 

November 1981 f' - ^ 

\ / 
FOR EARLY DOMESTIC DISSEMINATION 


iW\SA 

National Aeronautics and 
Space Administration 

Langley Research Center 

Hampton, Virginia 23665 


Because of its significant early commercial potential, this 
information, which has been developed under a U S Gov- 
ernment program, is being disseminated within the United 
States in advance of general publication This information 
may be duplicated and used by the recipient with the ex- 
press limitation that it not be published Release of this 
information to other domestic parties by the recipient shall 
be made subject to these limitations 

Foreign release may be made only with prior NASA ap- 
proval and appropriate export licenses This legend shall 
be marked on any reproduction of this information in whole 
or in part 

Review for general release November 30, 1983 



NF01426 



All Blank Pages 

/ 

Intentionally Left Blank 
To Keep Document Continuity 



TABLE OF CONTENTS 


SECTION PAGE 



NOMENCLATURE 

ii 

I 

INTRODUCTION 

1 

II 

THEORETICAL DEVELOPMENT 

3 


Droplet Dynamics 

3 


Evaporation Model 

5 


Wake Flowfield 

7 


Basic Wake Model 

8 


Tunnel Model 

9 

■ 

Propeller Model 

9 

\ 

Crosswind Model 

11 

- 

Distribution and Drift Model 

11 

Ill 

CODE VALIDATION 

15 

IV 

COMPUTER PROGRAM MANUAL 

18 


General Description 

18 


Input Data 

19 


Output 

23 


Error Messages 

25 


Sample Cases 

27 

V 

SUMMARY 

28 


REFERENCES 

29 


FIGURES 

30 


APPENDIX A - SAMPLE CASES 

37 


APPENDIX B - PROGRAM LISTING 

61 


i 




yfe- 7/^39 



NOMENCLATURE 


A 

b 

C 



F 

Fr 

g 

h 

K 


Mv 


P 

Pf 

Py 

PKl to PK4 

Q 

Qp 

r 

R 

Sc 

T 


t 

U 


u 

V 

^cw 

X, y, z 
Xo> yO)2o 


Aircraft aspect ratio 
Aircraft semispan, m 
Concentration, volume/semispan'^ 

Drag coefficient, D/qs 
Lift coefficient, L/qs 
Diffusivity, cm2/sec 
Propeller thrust, N 
Froude number, U//bg 
Gravitational constant, m/sec^ 

Aircraft altitude, semispans 
Inertia parameter, oS^u/lSbp 

Mean molecular weight of gas mixture in transfer path 
Mean molecular weight of evaporating material 
Mass, Kg 

Pressure, Kg/m^ (Ib/in^) 

Partial pressure of air, N/m2 
Vapor pressure of air, N/m2 
Propeller constants 

Nozzle flow rate, volume/ semispan along flight path 

Propeller torque, N-m 

Radial coordinate, in propeller radii 

Droplet Reynolds number p6U|u - n | / y 

Droplet free stream Reynolds number, p6U/y 

Propeller radius, m 

Schmidt number, pDy/y 

Temperature, K(°R) 

Time, secs 

Free stream velocity, m/sec 

Local flowfield velocity, dimensionless with U 
Cummulative volume percent 
Crosswind velocity component, m/sec 
Cartesian coordinates 

Initial droplet position (nozzle location) , semispans 


a 

f^s 

r 

6 

‘^m 

n 

y 

p 

T 


Droplet density, Kg/m^ 

Standard deviation, microns 

Aircraft circulation strength, m^/sec 

Droplet diameter, microns 

Volume median droplet diameter, microns 

Droplet position, dimensionless with b 

Air absolute viscosity, Kg/m-sec 

Air density, Kg/m^ 

Dimensionless time, Ut/b 


Superscripts 

Derivative with respect to t 
V ector quantity 

Subscripts 

g Final value at ground plane 

wb Wet bulb condition 

db Dry bulb condition 
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I . INTRODUCTION 


The use of aerospace technology to improve aerial applica- 
tions promises to yield significant improvements in the effi-^ 
ciency of operation, the environmental acceptability, and many 
other areas. One area of study is the interaction of the spray 
droplet with the aircraft wake. This interaction improves the 
overall efficiency of the operation by increasing the swath 
width, but also results in the drift of particles suspended in 
the wake to off target sites. The propeller slipstream and other 
nonuniform portions the wake also lead to nonuniform deposits of 
material inside the swath. This research uses an analytical 

t. 

approach to study this wake-particle interaction. Ultimately 
the goal of such a research effort will be to tailor the aircraft 
wake and dispersal system to produce a wide uniform distribution 
of material with minimum drift. 

Reed^ in 1954 presented a two-dimensional analysis of aerial 
applications. This report, aided by the techniques for trajec- 
tory analysis developed in the aircraft icing program, established 
the fundamentals of the method. This work was hampered by the 
limited computational capabilities of the time and Reed presents 
only a limited number of trajectories for his simple model. 

n 

Bragg improved this model by performing a three dimensional 
analysis including propeller slipstream effects Trayford and 
Welch^ have made similar improvements including the effect of 
crosswind on the distribution. The present method is based on 
reference 2 but includes many improvements which have been 
developed by the author since the publication of reference 2. 



The goal of this research was to document a method, and 

develop a computer code to simulate the dispersal of liquid 

sprays from agricultural aircraft. The program was to require 

small computation times so mathematical models and numerical 

procedures were kept as simple as possible while still modelling 

the fundamental physical phenomena. The program was written in 

such a way that a user could modify the v?ake model easily, v/hile 

still using the particle dynamics model and the distribution and 

drift numerical procedures. The user may provide a flowfield 

code which can be inserted into this program by use of the USERV 

subroutine. This allows for the analysis of more complicated 

flowfields which. may. contain fuselage vortices, vortices off of 

/ 

flaps due to variable span loading, winglets, or more complicated 
flowfields due to unconventional aircraft configurations 

The computer code consists of three main parts; the droplet 
dynamics and evaporation models, the aircraft wake models, and 
the distribution and drift prediction method. The theoretical 
basis for these models is derived in Section II and the method 
is compared to other results in Section III. Section IV is the 
computer program users' manual containing a description of the 
code and the input, output, and error messages. In the Appendices 
are sample cases and a complete listing of the code. 
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II. THEORETICAL DEVELOPMENT 


In this section the mathematical models used for the droplet 
dynamics and the v?ake flowfield mil be described. The numerical 
procedures used to solve these equation vzill also be presented. 

Droplet Dynamics 

By application of Newton's Second Law of Motion to a water 
droplet the differential equation describing the particle tra- 
jectory may be derived. Here particles in the 100 to 500 micron 
size range are considered. These particles are nearly spherical 
in shape‘s permitting the use of standard sphere drag data. In 
addition the analysis does not include the effect on the trajec- 
tory of atmospheric and aircraft induced turbulence, liquid drop- 
let deformation and internal circulation, magnus forces, multiple 
particle interaction, and electric charge. For the aerial appli- 
cation problem these forces are considered negligible. 

The equation of motion for a non-evaporating single particle 
moving in the wake of an aircraft can then be written as 

d^X _ _ _ 

m(^) = D + P + Ma + B + mg (1) 

The apparent mass, Ma, the pressure gradient term, P, and the 
Bassett force, B, become important only if the density of the 
particle is of lower, or similar order of magnitude as that of 
air^ . In addition, the apparent mass and Bassett force terms 
could also become significant if the particle experiences a large 
acceleration. The particle could experience a sizeable accelera- 
tion when initially injected into the flow. This acceleration 
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is, however, of short duration and in the streamwise direction 
so as to have little effect on the final spanwise location of 
the particle. Therefore, noting that the density of water is 
much greater than that of air, and ignoring the initial accelera- 
tion effects; the apparent mass, the pressure gradient, and the 
Bassett force terms may be dropped from equation (1) . 


Now writing this equation in nondimensional form, 
•• Cdr^_ 1 ^ 


Kn = 


24 


(u - n) + 


( 2 ) 


where the nondimensional numbers K, the inertia parameter, and 
F^, the Froude number are 


K = 


ISbp 



(3)a,b 


An additional nondimensional number, Ry, the free stream droplet 
Reynolds mamber 




p6U 


(4) 


CqR 

results from the term in equation (2). The droplets studied 

here experience relatively low Rejmolds numbers well within the 
range of zero to one thousand. Several curve fits to the standard 
sphere drag curve are available in this range and all yield similar 
results. Here the equation by Langmuir^ 

^ = 1. + 0.19R°-^3 + 2.6 X 10->R^-^^ (5) 


is used to calculate the sphere drag when computing the particle 
trajectories. Note that here R is the Re5molds number experienced 
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by the droplet based on the particle velocity relative to the 
surrounding fluid 



The trajectory of a particle is then governed by equation (2) 
which is a non-lfnear, second order, ordinary differential equa- 
tion which must in general be solved numerically The equation 
may be reduced to first order and separated into its three com- 
ponents to form a system of six simultaneous differential equa- 
tions Given the six initial conditions of particle position 
and velocity when first injected into the aircraft flowfield, 
the trajectory may be calculated using a step integration method. 
The program uses a variable step size, predictor-corrector 
scheme suitable for stiff systems^' T^en compared against 
non-stiff methods this scheme shows a considerable decrease in 
computational time to achieve the same accuracy. The trajectory 
is terminated x^hen the particle intersects the ground plane or 
becomes entrained in a vortex. 

Evaporation Model 

For small liquid droplets under roughly 150 microns in dia- 
meter evaporation effects often become significant As a result 
of evaporation the ground distribution of material is reduced 
and the material lost due to drift is greatly increased. To 
include the effect of evaporation, the analysis of Goering^ is 
used. 

Since the mass of the particle is no longer constant the 
differential equation (1) must be modified. Having already 
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incorporated the assumptions described early it becomes 


/ d^x 
mC- 


dt^ 


) + 




dm 

dt 


= D + mg 


( 7 ) 


Nondimensionalizing, this equation becomes 


K(Sii) - 

dx2 




^i) 


3k,d(6/b).^ 1 

(T/h) ^ dx ^ p^2 



O 


( 8 ) 


Equation (8) is identical to equation (2) except for the addition 
of the second term on the right hand side which is due to the ^ 


term in equation (7). An additional variable, 6, the particle 
diameter is introduced in equation (8) and another equation must 
be added to the system to allow a solution. 

The expression for the rate of change of the particle dia- 

9 

meter with respect to time is 


>d(5/b), _ ,2, o l/3„l/2, 


+ 0 6 Se- 


R 


(9) 


The terms on the right hand side may be evaluated in the follow- 
ing way. 

It is assumed that the water droplet is diffusing into the 
air, therefore, the molecular weights are = 29 0 and My = 18.0. 
The diffusion coefficient for xcater vapor into air is given as a 

9 

function of temperature 

Dy = 5.28 X 10“^ (10) 


The term AP/Pf is shown by Goering^ to be 


AP 

P^ 


l*swb 

P 

^atm 


- P 


( 11 ) 
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Brooker^^ has developed a mathematical model for the psychro- 
metric chart which is useful in evaluating equation (11). Input 
to the computer code are the atmospheric pressure, and 

the wet and dry bulb temperatures, and The saturation 

vapor pressure, at wet bulb temperature, T^|^, is given by 

■;/, 12301.688 

Pswb = _ 5.16923 In(T^b) (12) 

and the vapor pressure, Py, is 


V ^swb 0.62194(1075.8965 - 0.56983(T^b " ^91.69)) 


+ 


0 . 2405 (Pswb “ Patm) 


CPdb'lwb) 

(13) 


In equation (12) and (13) the temperatures are in °R and the 
pressures are in psi. Equation (12) is valid only for temperatures 
above 32°F. Note that by definition RH = Py/Ps > H ^1^^ relative 
hunidity, RH, is available when is not. 

The particle trajectory including the evaporation effect is 
calculated in the same manner as before, here using the differential 
equation (8) and adding to the system equation (9) . Additional 
complications due to droplet evaporation which arise in the dis- 
tribution and drift calculations will be discussed latter 


Wake Flowfield 

The flowfield model incorporated in the computer code is a 
horseshoe vortex system in ground effect. This simple model in 
either its 2-D or 3-D forms provides the wake velocities with a 
minimum of computer time, allowing for the small execution times 
characteristic of this program. In addition, a simple propeller, 
crosswind, and tunnel flowfield models are included in the code. 
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If these simple models are inadequate, the user may supply an 
alternate flowfield model by providing the proper subroutine. 
Here the flowfield models included in the code are discussed. 


Basic Wake Model : The wake of an aircraft in ground effect 

may be modelled in three dimensions using a horse shoe vortex 
system. This model is reasonably accurate provided that, the 
wake roll-up occurs in a time much less than that of the particle 
trajectory, the particle trajectory intersects the ground before 
this system decays appreciably, and the particle remain outside 
of the viscous vortex cores and in the potential flow region. A 
two dimensional model can be formed consisting of only two doubly 
infinite straight trailing vortices and their reflection through 
the ground plane, figure 1. The strength of these vortices is 

calculated assiaming an elliptical loading 

4CLUb 


r = 


7TA 


(14) 


The motion of these vortices is given by Bragg and the velocity 
can be calculated using the Biot-Savart law^. This model can be 
extended to three dimensions by adding a bound vortex of strength 
r on the aircraft quarter chord line, figure 2. The Biot-Savart 
Law is used to calculate the velocity induced in the x direction^ 
and the y and z velocities are the same as in the two dimensional 
model . 


It should be noted that the two dimensional model uses Reed's^ 
correction for the bound vortex. If the user inputs zero for the 
initial droplet velocity, the program assigns to the droplet an 
initial velocity equal to the velocity the particle would have if 
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falling in free air in a constant flowfield with the velocity 
components equal to that existing in the 2-D model at the point 
of particle injection. This correction is ignored if instead 
non-zero initial droplet velocities are input to the computer 
program. 

Tunnel Model - The two dimensional model of reference 2 has 
been included in this program. This permits the calculation of 
droplet trajectories from an aircraft model in a wind tunnel to 
evaluate the influence of the tunnel walls and the floor. The 
tunnel ceiling is assumed to have a negligible effect on the 
droplet trajectories since a model in ground effect will be in 
close proximity to the tunnel floor. The simple geometry of this 
model is shown in figure 3. A confonnal mapping method is used 
to determine the velocity at a point in the tunnel due to the 
two dimensional vortex system. The complex function describing 
the velocity potential and stream function in free air is trans- 
formed to the plane with infinite reflections of this vortex 
system in the y direction by using a complex sine mapping func- 
tion. The vortex position is determined by step int-egration using 
a truncated series to exnress the induced velocities on a vortex 
filament. The derivative of this complex function with respect 
to the complex particle position, yields the y and z velocities 
at any point in the tunnel. 

Propeller Model : A propeller slipstream is characterized 

by a rotation in the y-z plane due to the torque and an increased 
axial velocity due to the generated thrust. Functional forms for 
the axial and rotational velocities as a function of radial 
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distance from the propeller hub have been assumed. These equa- 
tions are of a form compatible with existing data on the thrust 
and torque loadings on general aviation propeller blades. The 
equation for the axial velocity is 

Vaxial = PK1(1. - r) + PK2(r) sin(rr) (15) 

and the rotational velocity is given by 

Vrot = PK3(1. - r) + PK4(r) sinCrr) (16) 

Figure 4 shows the velocity distributions for a typical propeller 
model. PKl and PK3 govern the velocities at the centerline while 
PK2 and PK4 control primarily the maximum velocities. For most 
applications PK3 = 0 and PKl is set equal to some small fraction 
of PK2. Using equations (15) and (16) the slipstream can be 
integrated to determine the thrust, F, and the torque, Qp. The 
thrust is 

F = (RP)2ttp [2U(0.166667PK1 + 0.189304PK2) 

+ 0. 0833333(PK1)2 + 0. 129006(PK1) (PK2) (17) 

+ 0. 0870040 (PK2)2j 

and for the torque, 

= 2ttp (RP)^ [ 0 . 033333 (PKl) (PK3) + 0. 036657 [(PKl) (PK4) 

+ (PK2)(PK3)] + 0. 0507039 (PK2) (PK4) + 0 083333U(PK3) 

+ 0.124801U(PK4)] (18) 

Fquations (17) and (18) are then used to determine the constants 
PKl, PK2, PK3 , and PK4 knowing the propeller thrust and torque 
This method assumes that the propeller slipstream extends dox-m- 
stream from the propeller disc plane to infinity x-^ith no dissipation. 
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In reality the rotational velocities are noticeably reduced 
due to viscous dissipation and the straighting effects of the 
main wing, fuselage, landing gear, horizontal, and vertical 
stabilizers. This reduction in rotational velocity is strongly 
dependent on the particular aircraft being modelled. It is 
recommended therefore that Vrot> that is PK3 and Pk4, be reduced 
from those calculated from equation (16) . Although little experi- 
mental data is available on which to base this, a reduction to 
75 percent of the calculated values seems reasonable from pre- 
liminary studies of conventional low wing monoplanes. 

Crosswind Model : The component of the wind perpendicular 

to the flight path and parallel to the ground plane is modelled 
by the equation^^ 

Vcw = VCW(z)0-25 (19) 

Here VCW is the crosswind component in m/sec 1 meter above the 
ground and z is the height above the ground plane in meters. The 
crosswind not only contributes a velocity in the y direction, but 
also causes a lateral translation of the entire aircraft wake 
system. This is modelled in the program by including the first 
order differential equations for the positions of the trailing 
vortices and propeller slipstream with the system of equations 
for the droplet trajectory and solving this enlarged system by 
the same step integration technique. 

Distribution and Drift Model 

By combining several droplet trajectories from a single nozzle 
location the ground deposition from the nozzle can be calculated 
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and the drift estimated. This procedure is straight forward 
for the case of no crosswind, evaporation, or propeller and was 
first developed for the aerial spray problem by Reed^ . Here an 
improved form of Reed's analysis is used for the simple case and 
major changes are made to handle the more difficult case of the 
unsymmetric and often double valued distributions. 

The concentration, C, from a single nozzle can be \i7ritten 
as 

= " S % ( 20 ) 

By superposition the concentration from any number of nozzles 
can be found The derivative d6/dyg can be found from the plot 
of the particle diameter as a function of the final lateral 
position where the particle intersects the ground plane. These 
points are cubic splined to obtain d6/dyg. This value must then 
be multiplied by dV/d6, a function of 5 and the particular nozzle 
used, which represents the volume of material released by the 
nozzle as a function of particle diameter. Here a normal dis- 
tribution of material about some volume median diameter , 5ni, is 
assumed xjhich gives 


dV ^ - 

/2tt 0g 



( 21 ) 


where Og is the standard deviation The Q in equation (20) is 

the volume flow rate of the nozzle The concentration of material 

on the ground is then given by equation (20) and all that remains 
is to determine the drift. 
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The drift is estimated by an extrapolation procedure which 
finds the smallest particle, within some maximum allowable error, 
which reaches the ground. The amount of material in the distri- 
bution from particles smaller than this drifted particle is found 
by integration of the normal distribution function. The allow- 
able error is a function of the number of standard deviations the 
particle size is from 6^ to control the error in terms of volume 
of material. Particles which encircle the wing-tip vortices are 
considered as drift and their trajectory calculations are ter- > 
minated. In addition, off target distribution is also considered 
drift The user can select a lateral semispan distance beyond 
V7hich the material is considered to be drift. 

This procedure must be modified for unsymmetrical distri- 
butions or cases including droplet evaporation. As a result of 
the propeller slipstream effects or the crosswind model the cal- 
culation of d5/dyg can become more difficult. Here it is not 
uncommon for several particle sizes to intersect the ground at 
the same yg location, thus making the function double valued 
The double valued and often erratic behavior of the curve makes 
the use of the cubic spline very difficult. Many points would 
be required to obtain a reasonable fit VJhen these complications 
arise the program abandons the cubic spline and instead uses 
linear interpolation between the calculated trajectory points. 

In regions of large changes in slope the program automatically 
inserts additional particle trajectories using again an error 
narameter weighted from the volume median diameter. This method 
provides good results in a minimum of computer time. Here the 
drift is estimated in the same manner as before. 
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When the droplet evaporation is considered the mass lost 
through evaporation must be accounted for. This results in an 
increase in material lost due to drift and a decrease in the 
ground deposition. The final droplet diameter, 6£, is deter- 
mined as a function of the initial diameter, 6 q, by using a cubic 
spline on the results of the trajectory calculations. The dis- 
tribution is then reduced by the ratio of In addition 

to the material which does not reach the ground or lies outside 
the target area, the amount of material which evaporates is 
determined by integration and included in the total drift. By 
using this procedure accurate drift and distributions are deter- 
mined from the calculated particle trajectories. 
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III. CODE VALIDATION 


Initial validation of the computer code was conducted by- 
comparing the results of this method to the early analytical 
work of Reed^. Reed's method uses a two dimensional wake model 
without propeller, evaporation, or crosswind effects. The dis- 
tribution is calculated using a least-squares fit of a second 
order polynomial to the trajectory data. A comparison of the 
two methods is shown in figure 5. Here the circulation, 'i' , has 
been matched since the early work assumed a rectangular lift 
distribution. The comparison is quite good considering that 
modern digital computers were not available for Reed to make 
his trajectory calculations. 

Comparisons of the code to experimental data is difficult 
because of the number of variables involved which are difficult 
to control or measure. Recent small scale model tests of agri- 
cultural aircraft dispensing scaled particles have generated some 

19 13 

distribution data which may be used for comparison ’. Figure 
6 shows the measured distributions from three nozzles and the cor- 
responding calculated distributions. The model used was unpowered 
with a wing semispan of 0.914m moving at 20.6 m/sec at an altitude 
of 0.51 semispans. The particles used were glass microbeads with 

O 

a mean diameter of 125 microns and a density of 2.42 g/m . Since 
the test program was designed only to validate the scaling laws, 
no measurements of particle size distribution or nozzle volumetric 
flow rate was made. These parameters were estimated for this com- 
parison. Note also that the concentration reported in reference 12 
is in terms of the numbers of particles per unit area where this 
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method calculates concentration in terms of volume of material 
per unit area. This accounts for the difference in the areas 
under the curve and also for the larger measured concentrations 
to the right of the peak concentration where the particles are 
smaller in diameter. The difference in the distributions at the 
y-Q = 0.30 nozzle location is due in part to the fuselage vortex 
not accounted for in the present flowfield model. The comparison 
in figure 5 is however, quite good overall and supports the analysis. 

In figure 7, the oresent method is again compared to data from 
reference 12. Here the same model is used as in figure 5 but the 
altitude is 0.35 semispans. The lateral position of the nozzle is 
at 0.40 semispans and the particles have a mean diameter of 140 
microns and a density of 0.58 g/m^. These particles scaled to a 
200 micron water droplet for the full scale aircraft. Here again 
the comparison is quite good considering the differences in the 
definition of concentration. The amount of lateral transport of 
the material and the general shape and magnitude of the concentra- 
tion curves compare quite well. 

The present method was also compared qualitatively to the 
powered test conducted in reference 12. Again considering the 
influence of the fuselage vortex the comparisons were encouraging. 

The amount of lateral transport due to the propeller was approxi- 
mately the same if the calculated swirl velocity v/as reduced by 
25% as discussed earlier to account for straighting effects. 
Sufficient data were not available to make a detailed comparison 
of the distributions. 

Detailed comparisons of this method have been conducted 
with both analytical and experimental data. The method has 
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performed well over the range of parameters tested 
appears to do a good job of predicting the lateral 
droplets in the w^ake of aircraft in ground effect. 


The code 
transport of 



IV. COMPUTER PROGRAM MANUAL 


A computer program has been developed to predict the ground 
deposition resulting from the dispersal of water droplets into 
the wake of an aircraft in ground effect. A complete listing 
of this program has been provided in Appendix B. After a general 
description of the program, the input, output, error messages, 
and test cases will be presented in this section. 

General Description 

The computer code has been written in FORTPvAN and developed 
on the AARL Harris Computer System. The code is in single pre- 
cision and should execute on CDC systems with only minor changes. 
Conversion of the code to execute on IBM type machines should 
involve little more than a conversion to double precision and 
modifying the Hollerith codes. The code is approximately 2500 
card images in length and requires 15 to 45 seconds of execution 
time per single nozzle deposition on the Harris SLASH 6 Computer. 
The Harris /6 is slightly slower than a CDC Cyber 73 and about 
6 times slower than an IBM 370/168. Additional storage must be 
supplied to the program on octal file 10. 

A plot package has been provided which is completely contained 
in subroutine PLOTl. A single call to this subroutine occurs 
as the last executable statement in the main program. Since the 
plot subroutine used here may not be compatible to those avail- 
able on other systems, the program was written to allow for easy 
modification of the plot package. 

Certain limitations on the size of the analysis that can be 
conducted considering the present DO loop parameters and array 
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dimensions should be discussed. The maximum swath width that 
can be used is +10 semispans and the number of nozzles input can 
not exceed 100. The program is limited to 15 particle trajec- 
tories per nozzle including those determined internally to 
estimate the drift and improve accuracy in a region of a large 
change in lateral transport. Therefore no more than 10 particle 
diameters should ever be input and in general 5 to 7 is usually 
sufficient. No more than 1000 steps can be taken in the calcul- 
ation of a single particle trajectory. The number of particles 
calculated and the number of steps can be controlled by use of 
the user input error parameters. 

A flow chart including the major subroutines is shown in 
figure 8. 

Input Data 

The program reads the input from device 5 in the manner 
described belov7. Several options are available to the user 
with the acceptable combination of options sho™ in figure 9. 

In addition the user may supply the wake flowfield to be used 
by the inclusion of a subroutine USERV. This allows the pre- 
sent particle dynamics, including evaporation effects distribu- 
tion, and drift models to be used with a different or more 
sophisticated wake model. USERV is described by comment State- 
ments in the listing. It may use by inclusion of the proper 
COMMON statements or subroutine calls use other variables and 
sections of the code as desired. 

The inputs to the program are primarily physical variables 
in SI units. Some variables have been nondimensionalized where it 
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seemed the most logical approach It is hoped that this type 
of input will make the code more useful. It should be noted 
however, that the problem can be nondimensionalized as in 
reference 11. The number of independent parameters is much less 
than the number of physical variables input to this code. The 
inputs could have been in terms of completely nondimensional 
numbers such as Ry, K, and F^- as discussed earlier. The use of 
this code for sensitivity analysis can, therefore, be greatly 
simplified by considering the nondimensional parameters. 

Card 1- (Format 13A6) Title 

Card 2: (Format 15) 

NDROPS Number of trajectories to be run. 

Set = 1 if NDIST = 1 

Cards 3-8 repeated NDROPS times 

Card 3: (Format 912, 2X, E15.6) 

N2D = 0 For a 3D run 

1 For a 2D run 

2 For a 2D run using the user supplied 
USERV subroutine 

N3D = 0 For a 2D run 

1 For a 3D run 

2 For a 3D run using the user supplied 
USERV subroutine 

NEVAP =0 No droplet evaporation 

1 Droplet evaporation included 

NPR0P2 =0 No propeller effects 

1 Propeller slipstream model included 

NCW =0 No crosswind effects (set VCW = 0) 

1 Crosswind effects included 

NTUN = 0 No tunnel model 

1 Tunnel model included 

NDIST = 0 Single particle trajectory run 
1 Ground depositions calculated 
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NPRINT = 0 Short output option 

1 Long output option 

NPLOT =0 No plotting 

1 Plot distribution of trajectories 

EPS Step integration error control parameter, 

usually set equal to l.E-5. 

Card 4 included only if NEVAP = 1 

Card 4: (Format 3F10.6) 

O 

PA Atomspheric pressure, N/m 

TDB Dry bulb temperature, °C 

TUB Wet bulb temperature, °C 

Card 5 included only if NPR0P2=1 
Card 5; (Format 6F10.4) 

ZPR0P2 Z coordinate of propeller hub, semispans 

PKl to 4 Propeller slipstream constants as defined 
in Section II 

RP Slipstream radius (usually equal to pro- 

peller radius) , semispans 

Card 6 (Format 7F10.4) 

A Aircraft aspect ratio 

CL Aircraft lift coefficient 

B Aircraft semispan, m 

U Aircraft speed, m/sec 

2 

G Gravitational constant, m/sec 

ZYO Initial Y coordinate of trailing vortex, 

semispans 

ZZO Initial Z coordinate of trailing vortex, 

semispans 
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Card 7: (Format 2E20.6, 2F10.4) 

DA Air Density (ignored if NEVAP = 1), Kg/m^ 

VIS Absolute air viscosity (ignored if 

NEVAP = 1) , Kg/m-sec 

VCW Crosswind velocity (=0. if NCW = 0), m/sec 

D Tunnel width (ignored if NTUN = 0) , m 

Card 8 included only if NDIST = 0 

Card 8: (Format 8F10.5) 

DIA Particle diameter, microns 

DD Particle density, Kg/m^ 

X, Y, Z Initial particle position, semispans 

UD, VD, Initial particle velocity relative to 
WD the aircraft fixed coordinate system, 

nondimensional with respect to U 

Cards 9 to 16 included if NDIST = 1 

Card 9: (Format 615) 

NCOL Number of nozzles 

NROW Number of initial particle sizes input 

NDIAMN = 0 All nozzles the same DIAMN 

1 DIAMN input for each nozzle 

NSTDEV = 0 All nozzle have the same STDEV 

1 STDEV input for each nozzle 

NQ =0 All nozzles have the same Q 
1 Q input for each nozzle 

NZNOZ = 0 All nozzles have same Z coordinates 

1 Z input for each nozzle 

Card 10: (Format 8F10 5) 

DD Particle density, Kg/m^ 

XO Nozzle X coordinate (ignored if N2D = 1) , 

semispans 

UDO, VDO, Initial particle velocity relative to 
WDO the aircraft fixed coordinate system, 

nondimensional with respect to U 

SWIDTH One-half the swadth width calculated, 
semispans 
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DERR Constant by which default error limits 

are multiplied, normally set equal to one 

DWIDTH Material landing beyond + this value are 
considered drift, semispans 

Card 11: (Formant 8F10.5) 

YNOZ Y coordinate of each nozzle (WOOL values) , 

semispans 

Card 12. (Format 8F10.5) 

SIZE Particle sizes input to represent distri- 

bution (NROW Values) , microns 

Card 13: (Format 8F10.5) 

DIAMN Median particle diameter, microns 

Input one value if NDIAMN = 0 
Input NCOL values if NDIAMN = 1 

Card 14: (Format 8F10 5) 

STDEV Standard deviation of particle size 

distribution, microns 

Input one value if NSTDEV = 0 
Input NCOL values if NSTDEV = 1 

Card 15. (Format 8F10.5) 

Q Nozzle flow rate, unit volume/ semi span 

along flight path 

Input one value if NQ = 0 
Input NCOL values if NQ = 1 

Card 16 (Format 8F10 5) 


ZNOZ Z coordinate of nozzle, semispan 

Input one value if NZNOZ = 0 
Input NCOL values if NZNOZ = 1 

Cards 1-16 may be repeated to run several cases simultaneously. 
Output 

The program first outputs all the basic input data. The 
variables output are all in the same units as the corresponding 
inputs. After this initial listing, the output depends primarily 
on only the choice of NDIST and NPRINT. The program's calculated 
output is all dimensionless unless indicated otherwise. 
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For the NDIST = 1 case and NPRINT = 1, the distribution from 
each individual nozzle is listed. The output is self explanatory, 
note that concentration has units of volume per semispan squared. 
If more than one particle size lands at the same lateral position, 
YG, the diameter, DIA, column is omitted. The rest of the output 
is the same for both the NPRINT = 0 and 1 options A table 
including the input values for each nozzle is output along with 
the drift in percent of total for each nozzle. Next a trajectory 
summary listing all the trajectories calculated is included before 
the final distribution The final distribution is a compilation 
of the distributions from all nozzles and is output in 0.01 semi- 
span increments. The total drift given is in percent of total 
material and is weightea by the volume flow rate of each nozzle. 

As a check, the total material reaching the tar j et area is also 
given. This is determined by a simple trapezoidal integration of 
the final integration and is meant only as a check. The accuracy 
of this estimate will be low for complicated distributions. 

For the single particle trajectory calculations of NDIST = 0 
the output is simple for NPRINT = 0. Here the initial inputs are 
output as before, followed by the final value of the variables in 
their normal units. The final values include particle position, 
velocity, droplet diameter for NEVAP = 1, and the nondimensional 
time, T, required for the trajectory. This output is repeated 
for each particle trajectory. 

If the user selects NPRINT = 1, all the dependent variables 
and the independent variable, x, are output at each time step in 
the integration process. This can result in a great many lines 
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in some cases, Since so many combinations of variables are 
possible depending of the user options selected, no header is 
provided for this output. The variables are output in the fol- 
lowing manner. Column 1 is always the step number and 2 the 
nondimensional time, x. Columns 3-6 always contain, the lateral 
position, Y, the lateral velocity, VD, the vertical position, Z, 
and the vertical velocity, WD If the run is 3-D, columns 7 and 
8 contain the streamwise position, X, and the corresponding 
velocity, UD. If NEVAP = 1 the last column always contains the 
current droplet diameter. The columns between WD or UD and the 
droplet diameter contain the vortex and propeller slipstream 
position. If options are selected making the entire flowfield 
transport laterally, the columns will be Y and Z coordinate of 
the left, then right trailing vortex and then the Y and Z position 
of the propeller slipstream. If all these variables are not pre- 
sent, the initial values make it clear which one is being output. 

Error Messages 

The program provides several internal error messages and 
warnings. These messages and the most likely user actions 
needed to correct the problem are given here 


MESSAGE 

COMMENT 

Max Number of Runs Exceeded 
Before Drift Calculation Com- 
pleted (E) 

Trajectory Terminated, Propel- 
ler Slipstream Entrainment 

Number of Particle Traj ectories 
Exceeds 15, Reduce NROW or 
Increase DERR 

Particle Considered Drift 
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MESSAGE 


COMMENT 


DO Loop Parameter Exceeded in Number of Steps Exceeds 1000, 
SUB2D (E) Increase EPS 

Trajectory Completed, Z less Normal Trajectory Termination 

Than Zero 

Trajectory Terminated, Vortex Particle Considered Drift 

Entrainment 

DO Loop Parameter Exceeded in Number of Steps Exceeds 1000, 
SUB3D(E) Reduce EPS 

NVAR Incorrect ^E) Illegal Combination of Options 

Cubic Spline Problem Nozzle Cubic Spline Provides Poor Fit 

No DIA = (E) to Trajectory Data, Change 

Error Parameters or Input Par- 
ticle Diameters 

Insufficient Data To Calculate Not Enough Trajectories to 
Distribution (E) Cubic Spline Particle Diameters 

for NEVAP = 1 Case, Change 
Input Particle Diameter 

Insufficient Data To Estimate (Same as Above) 

Drift, Nozzle No. Drift Set 
Equal To Zero (E) 

WDR Exceeded DO Loop in INCON Error in 2D Initial Droplet 

(E) Velocity Calculation, Input Non 

Zero Initial Velocity, Often 
Due to Illegal Input 

Right Endpoint in INCON Incor- (Same as Above) 
rect (E) 

DO Loop Parameter in INCON (Same as Above) 

Exceeded (E) 


(E) ERROR (W) WARNING 
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Sample Cases 

Five sample cases have been provided in Appendix A. These 
cases are intended both as examples and as check cases after the 
program has been installed on a new system. The complete input, 
output, and plots are given for each of the cases. With the use 
of the input and output descriptions in this section, the cases 
should be self explanatory. 
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V SUMIIARY 


A method has been developed and o computer program written 
to predict the dispersal of aerial sprays from agricultural air- 
craft. The method has been compared to early analytical results 
and recent experimental studies and found to compare well with 
these results A user's manual for the program including sample 
cases and FORTRAN program listing have been provided The code 
was written to allow the flexibility to install a more sophisti- 
cated wake model to increase the usefullness of the program. 

The computer program provides a useful tool to aid in agricultural 
aircraft and dispersal system design as well as for aerial appli- 
cations research 
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Figure 1. Two Dimensional Wake Vortex Model 



Figure 2. Three Dimensional Wake Vortex Model 


30 



CONCENTRATION VOLITIE/UNIT AREA 



Lateral Distance, Semispans 


Figure 5. Present Ilethod Compared to that of Reed 
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Figure 6. Comparison of Theory to Experiment for Three Nozzle Locations 


CONCENTRATION 







PLOT 1 


PLOT 

SUBROUTINES 
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Figure 9. Program Option Available 















































































































APPENDIX A 


FIVE SAMPLE CASES 

All five sets of input data are listed first, then each 
sample output including the plot. 
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OSU / AARL 


AERIAL SPRAY PROGRAM 
last update 2/1/81 


ASP TEST Case i 


program control 


N?D= 0 

N3D= I 

NEVAP= 0 

NPROP2= 1 

NCW= 1 

NTUN= 0 

NOIST= 0 

NPRINT= 1 

NPLOT= 1 

EPS= O.lOOOOOE-OA 


PARAMETER VALUES 
A= 6.0000 

CL= O.GOOO 

B= 6.0960 

U= 53.0350 

G= 9. 8066 

ZYOs 

O.BOOO 

7Z0= 0.5000 

OAs 0.122600E+01 

VIS= 0.179320E-OA 

VCW= 2.5000 

n= 

0 .0000 

nn= loon. 0000 





PR0P2 INPUTS 








2PH0P2= 0.5000 

RP= 0.2000 

PK1 = 

1.7930 

PK?= 17.9300 

PK3 = 

O.QOOO 

PKA = 

6.0290 


INITIAL values 


x= 

0.000000 

Y = 

o.iooooo 

Z= 0.400000 

W0 = 

0.000000 

niA= 

200.0000 


1 

O.UOOlO 

O.loOOO 

0.00003 

0.40000 -0.00011 

0.00000 

2 

0.00020 

0. inooo 

0.00005 

0.40000 -0.00022 

0.00000 

3 

0.00116 

o.ioooo 

0.00030 

0.40000 -0.00124 

o.ooool 

4 

0.00212 

0 . 1 0000 

0.00054 

0. 40000 -0.00225 

0.00003 

5 

0.00308 

o.ioooo 

0.00077 

0.39999 -0.00323 

0.00007 

6 

0.00472 

o.ioooo 

O.OOl 16 

0.39999 -0.00487 

0.00015 

1 

0.00636 

o.lnool 

0.00153 

0.39998 -0,00644 

0.00027 

a 

O.OOflOl 

0. loool 

0.00189 

0.39997 -0.007:<7 

0.00043 

9 

0.00965 

o.loool 

0.00224 

0.39995 -0,00945 

0.00061 

10 

0.01185 

0. 1 0002 

0.00268 

0.39993 -0.01136 

0.00091 

11 

0,01405 

0.10002 

0.00310 

0.39990 -0.01321 

0.00126 

12 

0.01625 

0,10003 

0.00350 

0.39987 -0.01499 

0.00166 

13 

0,01845 

o.lnn04 

0.00388 

0.39984 -0.01671 

0.00211 

14 

0.02225 

0.10005 

0.00449 

0.39977 -0.01957 

0.00299 

15 

0.02605 

0.10007 

0.00505 

0.39969 -0.02231 

0.003/9 

16 

0.02986 

0.10009 

0.00557 

0,39960 -0.02494 

0.00512 

17 

0.03366 

O.looil 

0.00606 

0.39950 -0.02749 

0.00635 

la 

0.03746 

0. 10014 

0.00650 

0.39939 -0.02998 

0.00768 

19 

0.04210 

0.10017 

0.00701 

0.39924 -0.03294 

0.00944 

20 

0,04674 

0. 10020 

0.00747 

0.39908 -0.03504 

0.01134 

21 

0.05139 

0 , 10024 

0.00789 

0.39891 -0.03870 

0.01335 

22 

0.05603 

0. 10020 

0.00828 

0.39H72 -0,04152 

0,01548 

23 

0.06067 

0.10032 

0.00864 

0.J9H53 -0.04432 

0.01772 


UD= O.OOUOOO VO= 0.000000 


0.00142 

0.80001 

0.50000 -o.aoooo 

0.50000 

0.00001 

0.53000 

0.00283 

0.80002 

0.50000 -0,79999 

0.50000 

0.00001 

0.49999 

0.01625 

0.80009 

0.50000 -0.79994 

0.50000 

0.00007 

0.49996 

0.02937 

0.80016 

0.49999 -0,79990 

0.49999 

0.00013 

0.49993 

0.04219 

0 .80024 

0.49999 -0.79985 

0.49999 

0.00019 

0.49990 

0.06347 

0.80036 

0.49998 -0.79977 

0.49998 

0,00029 

0,49985 

0.08395 

0.80049 

0.49998 -0.79970 

0.49998 

0.00040 

0.49980 

0.10368 

0.80062 

0.49997 -0.79962 

0.49997 

0.00050 

0.49975 

0.12269 

0.80074 

0.49997 -0.79954 

0.49997 

0.00060 

0.49970 

0.14711 

0.80091 

0.49996 -0.79943 

0.49996 

0.00074 

0.49963 

0.17040 

0.80108 

0.49995 -0.79933 

0.49995 

0.00088 

0,49957 

0.19263 

0.80125 

0.49994 -0,79922 

0.49994 

0.00101 

0,49950 

0.21387 

0 .8014? 

0.49993 -0.79912 

0.49993 

0,00115 

0.A9943 

0.24838 

0.80171 

0.49992 -0.79894 

0.49992 

0.00139 

0.49931 

0.28043 

0.80200 

0.49991 -0.79876 

0.49991 

0.00162 

0.49920 

0.31025 

0,80229 

0.49989 -0,79858 

0.49989 

0,00186 

0.49930 

0.33808 

0.80259 

0.49988 -0.79839 

0.49988 

0.00210 

0.49896 

0.36411 

0.80288 

0.49987 -0.79821 

0,49987 

0.00233 

0.49884 

0.39372 

0.80324 

0.49985 -0,79799 

0.49985 

0.00262 

0,49870 

0.42117 

0.80359 

0.49983 -0.79777 

0.49983 

0.00291 

0.49856 

0 .44670 

0.80395 

0.49982 -0,79755 

0.49902 

0.00320 

0.49842 

0.47052 

0.80431 

0.49980 -0.79733 

0.49980 

0.00349 

0.49827 

0.49281 

0.80466 

0.49978 -0.79711 

0.49970 

0.00378 

0.49al3 
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o.ufafl/ 

0. 10043 

0.00945 

0.39794 -0.05157 

0.02406 

0.54516 

0.80550 

0.49974 -0,79652 

0.49974 

0,00454 

0.49775 

?fy 

0.07H97 

0. 10049 

0.00980 

0.39761 -0.05516 

0.02745 

0.55851 

0.80607 

0.49972 -0.79623 

0.49972 

0,00492 

0.49757 

27 

0.08506 

0.10055 

0.01011 

0.39727 -0.05873 

0.03099 

0.59031 

0.8065* 

0.49970 -0.79594 

0.49970 

0.00529 

0.49738 

2H 

0.09116 

0.10061 

0.01040 

0.39690 -0.06228 

0.03465 

0.61073 

0.80701 

0,49968 -0.79565 

0.49968 

0,00567 

0.49719 

29 

0.09726 

0. 10067 

0.01065 

0.39651 -0.06579 

0.03843 

0.62994 

0,80748 

0.49965 -0.79536 

0.49965 

0.00605 

0.49700 

30 

0.10709 

0. 10078 

0.01102 

0.39583 -0.07136 

0.04477 

O.558TO 

0.80823 

0.49962 -0.79489 

0.49962 

0.00666 

0.49670 

31 

0.11692 

0.10089 

0.01133 

0.39511 -0.07679 

0.05137 

0.68515 

0.80899 

0,49959 -0.79442 

0.49959 

0.00728 

0.49640 

32 

0.1267b 

u.inioo 

0.01159 

0.39433 -0.08203 

0.05823 

0.70963 

0,80974 

0.49955 -0.79396 

0.49955 

0,00789 

0.49610 

33 

0.13658 

0.10112 

0.01181 

0.39349 -0.08705 

0.06532 

0. 732-1 

0.81050 

0.*9952 -0.793*9 

0.49952 

0.00850 

0.49580 

34 

0. 1A641 

0.10123 

0.01199 

0.39261 -0.09179 

0.07262 

0.75378 

0.81125 

0.*99*B -0.79302 

0.499*8 

0.00911 

0.49550 

3b 

0.1S623 

0.10135 

0.01214 

0.39169 -0.09624 

0.08013 

0.77381 

0.81201 

0.*994S -0.79255 

0.49945 

0.00972 

0.49520 

36 

0.17090 

0.10153 

0.01231 

0.39023 -0.10229 

0.09169 

0.80153 

0.B1314 

0,*9939 -0.79185 

0,49939 

0.01063 

0.49475 

37 

0.18556 

0.10171 

0.01242 

0.38869 -0.10758 

0.10363 

0.82692 

0.01427 

0.*9934 -0.79115 

0.49934 

0.01154 

0.49430 

38 

0.20023 

0. 10190 

0.01249 

0.38708 -0.11213 

0.11593 

0.85024 

0.81539 

0,*9929 -0.790*5 

0.49929 

0.012*5 

0.49385 

39 

0.21489 

0.10708 

0.01251 

0.38541 -0.11596 

0.12856 

0.87167 

0.81652 
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0.54128 

0.00328 

0,00002 -0.01269 

14.40428 

0.99999 

1.91644 

0.46472 -0.15507 

0.46472 

0.78432 

0.22665 

138 

14.28627 

0.54129 

0.00321 

0.00000 -0.01268 

14.40594 

0.99999 

1.91657 

0.46472 -0.15500 

0.46472 

0.78440 

0.22664 

139 

14.28793 

0.54129 

0.00314 

-0.00002 -0.01268 

14.40760 

0.99999 

1.91670 

0.46472 -0.15492 

0.4647? 

0,78447 

0.22662 


avail TRAJECTORY COMELEtET). 7 LESS THAN ZERO •»•« 


FINAL VALUES 

X= H.AOSROI 
WD= -0.0126R4 


Y= n.S41290 

niA= Ron.nooo 


z= o.nooooo 
T^ 14.?n6?97 


UD= 0.999992 


VO= 0.003211 
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ASP TEST case 2 


PROGRAM control 

N2D= 1 n30= 0 

NTUN= 0 NDIsT= 1 


PARAMETER VALUES 

A= 6.0000 CL= 0.6000 

ZY0= O.aOOO 770= 0.5000 

0= o.ooon 00= 1000.0000 


OISTRinUTION INPUTS 

NCOL= 1 NR0w= 6 

NZNOZ= 1 00=1000.00000 

■P' wo= 0.00000 <;winTH= a.ooooo 


osu / aarl 


AERIAL SPRAY PROGRAM 
LAST UPDATE 2/1/81 


nevap= 0 
nprint= l 

NPH0P2= 0 
NPLOT= 1 

NCW= 

EPS* 

0 

O.IOOOOOE-OA 

8= 6. 0060 
DA= 0.122500E+01 

U= 53.0350 

VIS= 0.179320E-OA 

G = 

VCW* 

9.8066 

0.0000 

N0IAMN= 1 
X= 0.00000 

DFRR= 1.00000 

NSTOEY= 1 
UO:: 0.00000 

0WIOTH= 5.00000 

NQ= ] 
V0 = 

1 

0.00000 


) 



YNOZ=0,400 


YG 01^ concentration YG 

0.6AOOOOE + 00 0.4A10(iie + n3 0.110564E-03 0.650000E + 00 

0.670000E*00 0.3y97f,?E + 03 0.33455GE-02 0.680000E + 00 

0.700000E+00 0.3665rsE*03 0.309327E-01 0.710000E+00 

0.730000E+00 0.J39334E+03 O.I3S946E+00 0.740000E+00 

0.760000E+00 0.JI6843E*03 0 . 36798bt+00 0.770000E+00 

. 0.790000E*00 0.397l9'iE+03 0.719134E + 00 0.8000Q0E + 00 

' 0.820000E*00 0.380570E+03 0.112l6lt+01 0.830000E+00 

0.850000E + QO 0.<;66l80F + n3 0.l4gi«0t*0l 0.860000E + 00 

, 0.88UOOOE+00 0.253538E+03 0.177154E*01 0,890000E*00 

0.910000E + 00 0.':424?fiE*03 0.193848E + 01 0.920000E + 00 

0.940000E-*00 0 ,<;32b7fiE*03 0.I99723E + 01 0,950000E + 00 

, 0.970000E + 00 O.H237q< 3E + 03 0.198964t + 0l 0.980000E*00 

' O.lOOOOOE + 01 0.«;15941E*03 0.18817 «E*o1 0,101000E*01 

0.103000E + 01 0.<iU887?E + m 0.1757I&E + 01 0.104000E*01 

f 0.106000E + 01 0 ,2U24R<iE + n3 n.l6141flE + Ol O.lOOOOOE + 01 

0. 109000E + 01 0.196691E + 03 0.H6626E + 01 O.llOOOOE + 01 

0.112000E+01 0.1914n9E+03 0.1323l3E*01 0.113000E+01 

, 0.115000E+01 0.186572E+03 0.118905E+01 O.11600OE+01 

' 0.118000E+01 0.1821??E+03 O.l06597t+ol 0.119000E+01 

0.121000E+01 0.178013E+03 0.9544984+00 0.1Z2000E+01 

0.124000F+01 0.174204E+03 0.854457L+00 0.125000E+01 

' 0.127000E + 01 0.U06fi?E + O3 0.78522ZE + 00 0.128000E + 01 

0.130000E+01 0.1673R0E+03 0.685937E+00 0.131000E+01 

, 0.133000E+01 0.16427?E+03 0.fil5658E*00 0.134000E+01 

^ 0.136000E+01 0.161377E+03 0.553441E+00 0.137000E+01 

0.139000E+01 0.158657E+03 0.498379E+00 0.140000E+OI 

, 0.142000E+01 0.1560q7E+03 0.449639E+00 0.143000E+01 

4> ' 0.145000E+01 0.l53bRlE+03 0.408463E+00 0.146000E+01 

Ln 0.148000E+01 0.1bl399E+03 0.368174E+00 0.149000E+01 

0.151000E+01 0.14923PE+03 0.334066E+00 0.152000E+01 

‘ 0.154000E+01 0.147l9?E+n3 0.302624E+00 0.155000E+01 

0.157000E+01 0.14b2ft?F+03 0.273632E+00 0.158000E+01 

, 0.160000F+01 0.14144SE+03 0.247128E+00 0.161000E+01 

' 0.163000E+01 0.14173RE+03 0.223049E+00 0.164000E+01 

0.166000E+01 0.140138E+03 0.201271E+00 0.1fi7oOOE+01 

f 0.189000E+01 0.13Hb37E+03 0.181636E+00 0.170000E+01 

' 0.172000E+01 0.137234E+03 0.163969t+00 0.173000E+01 

0.175000E+01 0.13b9?4E+03 0.148094E+00 O.176O0OE+O1 

, 0.178000E+01 0.134701E+03 0.133837t+00 0.179000E+01 

' 0.181000E+01 0.133581E+03 0.121034E+00 0.182000E+01 

0.184000E+01 0.1324qOE+o3 0.109534E+00 O.ISSOOOE+Ol 

0.187000E+01 O.lJlbllE+03 0.991974E-01 0.188000E+01 

^ 0.190000E+01 0,130Sq?E+03 0.89R979E-01 0.191000E+01 

0.193000E+01 0.129738E+03 0.816006E-01 0.194000E+01 

0.196000F+UI 0.12894PE+03 0.744098E-01 0.197000E+01 

'■ 0.199000E + 01 0 . 1 Z 82 onE + 03 0.681652E-01 0.200000E + 01 

, 0.202000E+01 0.127503E+03 0.827181E-01 0.203000E+01 

0.205000E + 01 0.l2684qE + n3 (1 . 579458E-0 1 0.206000E + 01 

" 0.208000E+01 0.12b231E*03 0.537473E-01 0.209000E+01 

0.211000E+01 0.125ft47F+03 0.500387E-01 0.212000E+01 

0.214000E+01 0.12b0q3E+O3 0.46750UE-01 0.215000F+01 

0.217000E + 01 0.U4587F + 03 0.438225E-01 0.218000E + 01 

0.220000E+01 0.iZ4084E+03 0.4120684-01 0.221000E+01 

0.223000E+01 0.iZ35R4F+03 0.388614E-01 0.224000E+01 

0.226000E+01 0.123124E+03 0.3675104-01 0.227000E+01 

0.229000E+01 0.122682E+03 0.348454E-01 0.230000E+01 

0.232000E+01 0.1222574+03 0.331191E-01 0.233000E+01 

0.235000E+01 0.121848E+03 0.3l5b0lt-01 o .236000E+01 

n <3 + at • a a r- ft 1 


DIA concentration YG DIA CONCENTRATION 
0.426213E+03 0.408700E-03 0.660000E+00 0.412488E+03 0.126168E-02 
0.387928E+03 0.779777E-02 0.690000E+00 0.376896E+03 0.162809E-01 
0.356927E+03 0.b41859E-01 0.720000E+00 0.347861E+03 0.88482BE-01 
0.331299E+03 0.198073E+00 0.750000E+00 0.323714E+03 0.275513E+00 
0.309752E+03 0.474192E+00 0.780000E+00 0.303311F+03 0.592145E+00 
0.291380E+03 O.8bl925E+00 0.810000E+00 0.285844E+03 0.9B7170E+00 
0.275b40E+03 0.125221E+01 0.840000E+00 0.270738E+03 0.137&34E+01 
0.261762E+03 0.159689E+01 0.870000E+00 0.257562E+03 0.169039E+01 
0.249681E+03 0.184000E+01 0.900000E*00 0.245979E+03 0.189559E+01 
0.239012E+03 0.196920E+01 0.930000E+00 0.235731F+03 0.198850E+01 
0.229539E+03 0.199634 e+oi 0.960000E+00 0.226616E+03 0.198681E+01 
0.221084E+03 0.194560E+01 0.990000E+00 0.218467E+03 0.191&22E+01 
0.213502E+03 0.184330E+01 0.102000E+01 0.211147E+03 0.180153E+01 
0.206672E+03 0.171079E+01 0.105000E+01 0.204544E+03 0.166297E+01 
0.200492E+03 0.156483E+01 0.108000E+01 0.198561E+03 0.151536E+01 
0.194877E+03 0.141774E+01 O.lllOOOE+Ol 0.193117E+03 0.136998E+01 
0.189751E+03 0.127730E+01 0.114000E+01 0.188139E+03 0.123259E+01 
0.185048E+03 0.114675E+01 0.117000E+01 0.1B3566E+03 0.110572E+01 
0.180717E+03 0.l027b2E+01 0.120000E+01 0.179347E+03 0.99036SE+00 
0.176711E+03 0.9l990bE+00 0.123000E+01 0.175442E+03 0.886566E+00 
0.172995E+03 0.823551E+00 0.126000E+01 0.171315E+03 0.793816E+00 
0.169536E+03 0.737734E+00 0.129000E+01 0.168436E+03 0.711317E+00 
0.166308E+03 0.661557E+00 0.132000E+01 0.165279E+03 0.638143E+00 
0.163286E+03 0.594069E+00 0.135000E+01 0.162322E+03 0.573341E+00 
0.160452E+03 0.534334E+00 0.138000E+01 0.159546E+03 0.515991E+00 
O.lb7787E+03 0.481469E+00 0.141000E+01 0.1S6934f+03 0.465232E+00 
0.155276E+03 0.4J4664E+00 0.144000E+01 0.154471E+03 0.420280E+00 
0.152906E+03 0.393188E+00 0.147000E+01 0.152146E+03 O.3804J3E+OO 
0.150665E+03 0.356392E+00 O.ISOOOOE+Ol 0.149945E+03 0.345064E+00 
0.148543E+03 0.323321E+00 0.153000E+01 0.147861E+03 0.312838E+00 
0.146536E+03 0.292683E+00 0.156000E+01 0.145893E+03 0.283019E+00 
0.144644E+03 0.264523E+00 0.159000E+01 0.14403BE+03 0.2556B9E+00 
0.142864E+03 0.238838E+00 0.162000E+01 0.142295F+03 0.230813E+00 
0.141192E+03 0.215b42E+00 0.165000E+01 0.140659E+03 0.208284E+00 
0.139626E+03 0.194496E+00 0.168000E+01 0.139126E+03 0.187953E+00 
0.138159E+03 0.175537E+00 0.171000E+01 0.137691E+03 0.169650E+00 
0.1J6788E+03 0.158487E+00 0.174000E+01 0.136351E+03 0.153197E+00 
0.13b507E+03 0.143170F+00 0.177000E+01 0.135100E+03 0.13R419E+00 
0.134312E+03 0.12941SE+00 0.180000E+01 0.133932E+03 0.125150E+00 
0.1331994+03 0.117063E+00 0.1B3000E+01 0.132845E+03 0.113232E+00 
0.132162E+03 0.10696bE+00 0.186000E+01 0.131833E+03 0.l02b22E+00 
0.131197E+03 0.9598B2E-01 0.189000E+01 0. 1308914+03 0.92889TE-01 
0.130300E+03 0.870086E-01 0.192000E+01 0.130016F+03 0.842382E-01 
0.129466E+03 0.790880E-01 0.195000E+01 0.129201E+03 0.7669334-01 
0.1286894+03 0.7223114-01 0.198000E+01 0.1284424+03 0.7015144-01 
0.1279634+03 0.6626754-01 0.201000E+01 0.1277314+03 0.6445324-01 
0.1272814+03 0.6105774-01 0.204000E+0l 0.1270634+03 0.5946824-01 
0.1266394+03 0.5648704-01 0.2070004+01 0.1264334*03 0.5508864-01 
0.1260334+03 0.5246044-01 0.2100004*01 0.125839F+03 0.5122514-01 
0.1254604*03 0.4889904-01 0.213000E+01 0.1252754*03 0. 4780344-01 
0.1249154+03 0.4573664-01 0.2160004*01 0.1247394+03 0.4478144-01 
0.1243974*03 0.4291814-01 0.2190004+01 0.124229E+03 0.4204684-01 
0.1239024+03 0.4039694-01 0.2220004+01 0.123742F+03 0.3961554-01 
0.1234294+03 O.J81334E-01 0. 2250004+01 0.1232754+03 0.374303E-01 
0.1229754+03 0.360944 e-01 0.2280004*01 0.122828E+03 0.3545954-01 
0.1225394+03 0.3425124-01 0.231000E+01 0. 1223974+03 0.3367614-01 
0.1221194+03 0.3257964-01 0.234000E+01 0.1219834+03 0.3205694-01 
0.1217144*03 0.3105874-01 0.2370004+01 0.1215824+03 0.3058214-01 
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FINAL distribution 


YG 

concentration 

YG 

CONCENTRATION 

YG 

CONCENTRATION 

YG 

concentration 

YG 

concentration 

0.6A 

0.110S64E-03 

0.65 

0.408700E-03 

0.66 

0.126168E-02 

0,67 

0.334558E-02 

0.68 

0.779777E-02 

0.69 

0.1628U9E-01 

0.70 

0.309327E-01 

0.71 

0.541859E-01 

0.72 

0.884828E-01 

0.73 

0.135946E+00 

0.74 

0.l98U/TE+n0 

0.75 

0.275S13E*00 

0.76 

0.367965E+00 

0.77 

0.474192E+00 

0.78 

0,592145E+00 

0.79 

0,719134E+00 

0.80 

0.851925E+00 

0.81 

0.987170E4U0 

0.82 

0.112161E+01 

0.83 

0.125221E401 

0.64 

0.137634E+01 

0.85 

0.149180E401 

0.86 

0.159689E401 

0.87 

0.169039E*01 

0.88 

0.177154E+01 

0.89 

0.184UOo£+fll 

0.90 

0.189559E+01 

0.91 

0.193848E+01 

0.92 

0.196920E+01 

0.93 

O.19805OE+O1 

0.94 

0,199723E*0l 

0.95 

0. 199634E*01 

0.96 

0.198681E*01 

0.97 

0.196964E+01 

0.98 

0.194580E+01 

0.99 

0.191622E*01 

1.00 

0*188l78E+0l 

l.OI 

0.184330E+01 

1.02 

0.180153E+01 

1.03 

0.175716E+01 

1.04 

0.171U^9E+o1 

1.05 

0.166297E*01 

1.06 

O.161410E*O1 

1.07 

0.156483E401 

1.08 

0.151536E401 

1.09 

0.146626E+01 

1.10 

0.141774E*01 

1.11 

0.136998E*01 

1.12 

0.132313E+01 

1.13 

0.127730E+01 

1.14 

0.123259E+01 

1.15 

0.118905E+01 

1.16 

0.114675E*01 

1.17 

0.110572E+01 

1.18 

0.106597E+01 

1.19 

0.102752E+01 

1.20 

0.990365E+00 

1.21 

0.954498E*00 

1.22 

0.919905E+00 

1.23 

0.886566E+00 

1.24 

0.854457E+00 

1.25 

0.8235S1E+00 

1.26 

0.793816E*00 

1.27 

0.765P22E+00 

1.28 

0.737734E400 

1.29 

0.71U17E + n0 

1.30 

0.685937E+00 

1.31 

0.661557E+00 

1.32 

0.638143E+00 

1.33 

0.615658E+00 

1.34 

0.594U69E+00 

1.35 

0.573341E*00 

1.36 

0.553441E+00 

1.37 

0.534334E+00 

1.38 

0.515991E+00 

1.39 

0.498379E+00 

1.40 

0«481469E*00 

1.41 

0.465232E»00 

1.42 

0.449639E+00 

1.43 

0.434664E+00 

1.44 

0.420280E+00 

1.45 

0.406463E+00 

1.46 

0.39318BE+00 

1.47 

0.380433E+00 

1.48 

0.368174E+00 

1 .49 

0.356J92E+00 

1.50 

0.345064E+00 

1.51 

0.334066E+00 

1.52 

0.323321E+00 

1.53 

O.312038E+OO 

1.54 

0.302624E+00 

1.55 

0.292683E+00 

1.56 

0.283019E+00 

1.57 

0.273632E+00 

1.58 

0.264523E+00 

1.59 

0.255689E+00 

1.60 

0.247128E+00 

1.61 

0.238B38E+00 

1.62 

0.230813E+00 

1.63 

0.223049E+00 

1.64 

0.215542E+00 

1.65 

0.208284E400 

1.66 

0.201271E+00 

1.67 

0.194496E+00 

1.68 

O.107953E+OO 

1.69 

O.101&J6E+OO 

1.70 

0. 175537E+00 

1.71 

0.169650E+00 

1.72 

0.163969E+00 

1 .73 

O.I58487E+00 

1.74 

0.153197E+00 

1.75 

0.148094E*00 

1.76 

0.143170E+00 

1.77 

0.138419E+00 

1.78 

0.133837E+00 

1.79 

0.I29415E+0O 

1.80 

0.125150E+00 

1.81 

0.121034E*00 

1.82 

0.117063E+00 

1.83 

0.113232E+00 

1.84 

0.1095J4E+00 

1.85 

0. 105966E+00 

1.86 

0.102522E+00 

1.87 

0.991974E-01 

1.88 

0.959882E-01 

1.R9 

0.928897E-01 

1.90 

O.098979E-O1 

1.91 

0.870086E-01 

1.92 

0.842382E-01 

1.93 

0.816006E-01 

1.94 

0.790880E-01 

1.95 

0.766933E-01 

1.96 

0.744098E-01 

1.97 

0.722311E-01 

1.98 

0.701514E-01 

1.99 

0.681652E-01 

2.00 

0.662675E-01 

2.01 

0.644532E-01 

2.02 

0.627181E-01 

2.03 

0.610577E-01 

2.04 

0.594682E-01 

2.05 

O.579450E-O1 

2.06 

0.564870E-01 

2.07 

0.550886E-01 

2.08 

0.537473E-01 

2.09 

0.5246U4E-01 

2.10 

0.512251E-01 

2.11 

0.500387E-01 

2.12 

0.488990E-01 

2.13 

0.478034E-01 

2.14 

0.467300E-01 

2.15 

0.457366E-01 

2.16 

0.447614E-01 

2.17 

0.438225E-01 

2.18 

0.429181E-01 

2.19 

O.42O460E-nl 

2.20 

0.412068E-01 

2.21 

0.403969E-01 

2.22 

0.396155E-01 

2.23 

O.308614E-O1 

a. 24 

0.381334E-01 

2.25 

0.374303E-01 

2.26 

0.367510E-01 

2.27 

0.360944E-01 

2.28 

0.354595E-01 

2.29 

0.348454E-01 

2.30 

0.342512E-01 

2.31 

0.336761E-01 

2.32 

0.331191E-01 

2.33 

0.325796E-01 

2.34 

0.320569E-01 

2.35 

0.315501E-01 

2.36 

0.310587E-01 

2.37 

0.305821E-01 

2.30 

0.301195E-01 

2.39 

0.2967U5E-01 

2.40 

0.292346E-01 

2.41 

0.288111E-01 

2.42 

0.283996E-01 

2.43 

0.279996E-01 

2,44 

0.276107E_01 

2.45 

0.272325E-01 

2.46 

0.268645E-01 






TOTAL DRIFT= 0.5580J2E+01 


TOTAL= 0.RAA312E+02 
0AKGNOI STOP 




OSU / AARL 


aerial spray program 
last update 2/1/fll 


asp test case 1 
program control 


N20= 1 

N1D= 0 

nevap= 1 

NPROP2= 0 

NCWs 0 

NTUN= 0 

NOIST= 1 

NppiNTs 0 

NPLOT= 1 

EPSa O.IOOOOOE 

parameter values 

A= 6.0000 

CL= 0.6000 

6.0960 

U= 53.0350 

G= 9.8066 

ZY0= 0.6000 

0 = o.nooo 

Z70= 0.5000 

no= 1000.0000 

DA= 0.118901E+01 

VI5= 0.181330E-OA 

VCW= 0.0000 

EVAPORATION INPUTS 

PA= 1013S2.10 

TOR= 33.8R8900 

TWR= 18.333300 



DISTRIBUTION INPUTS 

NC0L= 1 

NROW= 6 

NDIAMN= 1 

NSTDEV= 1 

NQ= 1 

NZNOZ= 1 

00=1000.00000 

X= 0.00000 

UD= 0.00000 

VO= 0.00000 

WO= 0.00000 

SWlnTH= 4.00000 

OERR= 1.00000 

OwIOTH= 5.00000 


NOZZLE NO. TNOZ 

ZNOZ 

DIAHN 

STOEV Q 

DRIFT 

1 O.AUOO 

0.4000 

POO. 0000 

50.0000 1.0000 

21.2474 



SINGLE droplet LATERAL DISPLACEMENT 


NOZZLEl 

YG 

0.S66S34F+00 
0,7a6095E+00 
0.892304E*00 
0.108410E401 
0.155559E+01 
0.204485E+01 
0.235994E+01 
O.OOOOOOE^Ol 
O.OOOOOOE^Ol 
0 .OOOOOOE+01 
O.OOOOOOE+01 
O.OOOOOOE+Ol 
O.OOOOOOE401 
O.OOOOOOE+OI 
O.OOOOOOE+01 


1 

Ut A 

0.6«00(inE*03 
0.3UOOOOE+03 
0.2bOOonE*03 
o.<iuooonF*o3 
0. 1500()nE + 03 
0.1JlH<iAE*n3 
0.1273?QE+03 
o.uuooooE+ni 
O.UUOOOOE+01 
O.ODOOOOE+01 
O.OUOOOOE+Ol 
O.UUOOOOE+Ol 
o.ouoonnE+ni 
O.OOOOOOE+01 
o.ooooooE+ni 


OIAG 

0.S98786E+03 
0.297565E+03 
0.247031E+03 
0.I96103E+03 
0.Ia 3102E*03 
0.n93<*2E*03 
0. 109949E+03 
O.OOOOOOE+OI 

o.onnoooE+ol 

O.OOOOOUE+01 

O.OOOOOOE+01 

O.OOOOOOE+OI 

o,ooooooe+oi 

O.OOOOOOE+01 

O.OOOOUOE+01 



final distribution 


YQ 

concentration 

YG 

CONCENTRATION 

YG 

CONCENTRATION 

YG 

concentration 

YG 

concentration 

0.64 

0.124109E-03 

0.65 

0 

445961E-03 

0.66 

0.134425E-02 

0.67 

0.349339E-02 

0.68 

0.800428E-02 

0.69 

0.164709E-01 

0.70 

0 

309088E-01 

0.71 

0.535747E-01 

• 0.72 

0.866979E-01 

0.73 

0.132177E+00 

0.74 

0.1913U8E+00 

0.75 

0 

264595E+00 

0.76 

0.351662E400 

0.77 

0.451289E+00 

0.78 

O.S61524E+00 

0.79 

0.679842E+00 

0.80 

0 

803113E+00 

0.81 

0.928154E+00 

0.82 

0.105193E+01 

0.83 

0.117164E+01 

0.84 

0.128490E*01 

0.85 

0 

138974E+01 

0.86 

0.148466E401 

0.87 

0.156864E*01 

0.88 

0.164107E+01 

0.89 

0.170174E+01 

0.90 

0 

175078E+01 

0.91 

0. 178863E+01 

0.92 

0.181580E+01 

0.93 

0.183294E+01 

0.94 

0.1B4078E+01 

0.95 

0 

184014E+01 

0.96 

0.183185E4U1 

0.97 

0.181677E401 

0.98 

0.179573E401 

0.99 

0.I769blE+01 

1.00 

u 

173887E401 

1.01 

0.170452E+01 

1.02 

0.166711E401 

1.03 

0.162724E+01 

1.04 

0.1b8543E*nl 

1.05 

0 

1S4217E+01 

1.06 

0.149790E+01 

1.07 

0.145297E+01 

1.08 

0.140773E+01 

1.09 

0.136a44E*ol 

1.10 

0 

131730E+O1 

1.11 

0.127254E+01 

1.12 

0.122836E+01 

1.13 

0.I1B492E+01 

1.14 

0.1 14237E+01 

1.15 

0 

IIOOBIE+Ol 

1.16 

0.106032E+01 

1.17 

0.102097E+01 

1.18 

0.982809E+00 

1.19 

0.945863E+n0 

1.20 

0 

910151E+00 

1.21 

0.875679E+00 

1.22 

0.842445E+00 

1.23 

0.810440E+00 

1.24 

0.779645E+00 

1.25 

0 

750041E+00 

1.26 

0.721601E+00 

1.27 

0.694295E+00 

1.28 

O.660O93E+OO 

1.29 

0.642961E+00 

1.30 

0 

618865E+00 

1.31 

0.595770E+00 

1.32 

0.573638E+00 

1.33 

0.552436E+00 

1.34 

0.532126E+00 

1.35 

0 

512674E+00 

1.36 

0.494045E+00 

1.37 

0.476205E+00 

1.38 

0.459121E+00 

1.39 

0.44a760E+00 

1.40 

0 

427091E+00 

1.41 

0.412084E+00 

1.42 

0.397710E+00 

1.43 

0.383941E+00 

1.44 

0.370749E+00 

1.45 

0 

358108E+00 

1.46 

0.345994E+00 

1.47 

0.334383E+00 

1.48 

0.323252E+00 

1 .49 

0.312579E+00 

1.50 

0 

302342E+00 

1.51 

0.292523E+00 

1.52 

0.283103E+00 

1.53 

0.274062E+00 

1.54 

0.2653BAE+00 

1.55 

0 

257052E+00 

1.56 

0.249039E+00 

1.57 

0.241244E+00 

1.58 

0.233645E4U0 

1.59 

0.226239E+O0 

1.60 

0 

219023E+00 

1.61 

0.211994E+00 

1.62 

0.205151E+00 

1.63 

0.198490E+00 

1.64 

0.192008E+00 

1.65 

0 

185704E+00 

1.66 

0.179573E+U0 

1.67 

0.173614E+00 

1.68 

0.167824E+00 

1.69 

0.162auoE*00 

1.70 

0 

156739E+00 

1.71 

0.151439E+00 

1.72 

0.146296E+00 

1.73 

0.141308E+00 

1.74 

0.136471E+00 

1.75 

0 

131783E*00 

1.76 

0.127240E+00 

1.77 

0.122840E400 

1.78 

0.118579E+00 

1.79 

0.114455E*n0 

1 .80 

0 

110463E+00 

1.81 

0.106601E+00 

1.82 

0.102866E+00 

1.83 

0.992536E-01 

1.84 

0.9b7616E-01 

1.85 

0 

923869E-01 

1.86 

0.8912b9E-01 

1.87 

0.859754E-01 

1.88 

0.829325E-01 

1.89 

0.7999J9E-01 

1.90 

0 

771566E-01 

1.91 

0.7441 76E-01 

1.92 

0.717737E-01 

1.93 

0.692221E-01 

1.94 

0.667b99E-nl 

1.95 

0 

643842E-01 

1.96 

0.620922E-01 

1.97 

0.598811E-01 

1.98 

0.577482E-01 

1.99 

0.5569U9E-01 

2.00 

0 

537067E-01 

2.01 

0.517930E-01 

2.02 

0.499474E-01 

2.03 

0.481675E-01 

2.04 

0.464509E-01 

2.05 

0 

447997E-01 

2.06 

0.432341E-01 

2.07 

0.417527E-01 

2.08 

0.403504E-01 

2.09 

0.390a27E-ni 

2.10 

0 

377650E-01 

2.11 

0.365733E-01 

2.12 

0.354437E-01 

2.13 

0.3A3724E-01 

2.14 

0.333562E-01 

2.15 

0 

323918E-01 

2.16 

0.314761E-01 

2.17 

0.306064E-01 

2.18 

0.297799E-01 

2.19 

0.289942E-01 

2.20 

0 

282470E-01 

2.21 

0.275361E-01 

2.22 

0.268593E-01 

2.23 

0.262148E-01 

2.24 

0.256007E-01 

2.25 

0 

250153E-01 

2.26 

0.244571E-01 

2.27 

0.239245E-01 

2.28 

0.234160E-01 

2.29 

0.2293U4E-01 

2.30 

0 

224664E-01 

2.31 

0.220229E-01 

2.32 

0.215986E-01 

2.33 

0.211926E-01 

a. 34 

0.208038E-01 

2.35 

0 

204314E-01 








TOTAL DRIFT= O.PiF'tTAE + oa 

TOTAL= 0.a60ABlE+Od 
BAKGND: STOP 
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OSU / AARL 


aerial spray program 

LAST UPDATE 2/1/Bl 


ASP TEST case a 
PROGRAM control 


N20= 1 

N3D= 0 

NFVAps 0 

NPROP2= 0 

NCW= 0 

ntun= 1 

nDist= I 

NPRINTa 0 

NPLOT= 1 

EPS* O.IOOOOOE' 

parameter values 

A= 6.0000 

CL= 0.6000 

a= 6.0960 

U= 53.0350 

G= 9.8066 

ZYO= 0.8000 

D= 36.S760 

7Z0= O.SOOO 

nn= 1000.0000 

Oa= 0.122S00E+0I 

VIS= 0.179320E-0A 

VCHs 0.0000 

DISTRIBUTION INPUTS 

NCOL= I 

NR0R= 6 

NDIAMN= 1 

NSTD£V= 1 

NQ= 1 

NZNOZ= 1 
W0= 0.00000 

00=1000.00000 
SWlPTHs A.OOOOO 

X= 0.00000 

DEPR= 1.00000 

UD= 0.00000 

0WIDTH= 5.00000 

VD= 0.00000 

NOZZLE NO. YNOZ 

znoz 

DIAMN 

STDEV Q 

DRIFT 

1 O.AUOO 

O.AOOO 

200.0000 

50.0000 1.0000 

S.5777 


single droplet lateral ntSPlACEMENT 


NOZZLE! 1 

YG 

0.567303E+00 

0.785011F+00 

0.8a8<JO2E*00 

0.107310E+01 

O.ISOOOlF+Ol 

0.I91081F+01 

0.?A8190E*01 

O.OOOOOOE+01 

O.OOOOOOE-^Ol 

O.OOOOOOE+01 

O.OOOOOOE'*’01 

o.oaaoooE-fOi 

O.OOOOOOE-fOl 

O.OOOOOOF+01 

O.OOOOOOE+01 


UIA 

0 •*>uoonnE*03 
o.30oonnE+n3 
O,25n0onE*O3 
o,2ooonnE*03 
o,i5oonnE+03 
0. 13029PE+03 
0. I20A17E+O3 

o.uuoonoE*oi 

O.UU0OonE*Ol 

OiOuoonoE+oi 

O.OOOOflOE+Ol 

o.ouoonnE+01 

o.uuooonF+01 

O.UOOOOOE*01 

o.ouooooE+ni 



FINAL DISTRIBUTION 


YG 

CONCtNTHATinN 

YG 

CONCENTRATION 

YG 

concentration 

YG 

concentration 

YG 

concentration 

0.64 

0.117847E.03 

0.65 

0 

•831420E-03 

0.66 

0.132101E-02 

0.67 

0.347890E-02 

0.68 

0.806143E-02 

0.69 

0.167483F__o1 

0.70 

0 

.316871E-01 

0.71 

0.553080E-01 

0.72 

0.900375E-01 

0.73 

0.137970E+00 

0.74 

0.200565t*n0 

0.75 

0 

.27843SE+00 

0.76 

0.371242E+00 

0.77 

0.477720E+00 

0.78 

0.595B02E+00 

0.79 

0.722T83E+OO 

0.80 

0 

.855J96E*00 

0.81 

0.990285E+00 

0.82 

0.112419E+01 

0.83 

0.125410E+01 

0.84 

0.137738E+01 

0.85 

0 

.189188E+01 

0.86 

0.1S9591E+01 

0,87 

0.168829E+01 

0.88 

0.176829E+01 

0.89 

0.183559E*nl 

0.90 

0 

.189011E*0l 

0.91 

0.193209E+01 

0.02 

0.19620BE+01 

0.93 

0.1980B3E+01 

0.94 

0.19(!9l8E + nl 

0.95 

0 

.1988096*01 

0.96 

0.197851E+01 

0.97 

0.196144E+01 

0.98 

0.193781E+01 

0.99 

0.190854E*01 

1.00 

0 

•1874S1E+01 

1.01 

0.183650E+01 

1.02 

0.179526E+01 

1.03 

0.175145E+01 

1.04 

0 . 1 70567E*nl 

1.05 

0 

.165845E+01 

1.06 

0. 161025E+01 

1.07 

0.156150E+01 

1.08 

0.151260E+01 

1.09 

0.146A03E+01 

1.10 

0 

•141601E*01 

1.11 

0.136872E+01 

1.12 

0.132230E+01 

1.13 

0.127688E+01 

1.14 

0.123253E+01 

1.15 

0 

.118933E*01 

1.16 

0.114733E+01 

1.17 

0.110657E+01 

1.18 

0.106707E+01 

1.19 

0.102884E+01 

1.20 

0 

.991886E*00 

1.21 

0.956193E+00 

1.22 

0.921753E+00 

1.23 

0,888547E*00 

1.24 

0.856552E+00 

1.25 

0 

.825743E+00 

1.26 

0.796091E+00 

1.27 

0.767564E+00 

1.28 

0.740131E+00 

1.29 

0.713757F+00 

1.30 

0 

.688409E+00 

1.31 

0.664052E+00 

1.32 

0.640652E+00 

1.33 

0.618174E+00 

1. J4 

0.596bO4E+n0 

1.35 

0 

.575848E*00 

1.36 

0.555935E+00 

1.37 

0.536811E+00 

1.38 

0.518445E+00 

1.39 

0.500807E+00 

1 .40 

0 

.4H3868E+00 

1.41 

0.467598E+00 

1.42 

0.451970E+00 

1.43 

0.436958E+00 

1.44 

0.422535E+nO 

1.A5 

0 

.408677E+00 

1.46 

0.395361E+00 

1.47 

0.382563E+00 

1.48 

0.370261E+00 

1.49 

0.358435E+00 

1.50 

n 

.347063E+00 

1.51 

0.336030E+00 

1.52 

0.325237E+00 

1.53 

0.314697E+00 

1.54 

0.304418E+00 

1 .55 

0 

.294405E+00 

1.56 

0.284664E+00 

1.57 

0.275197E+00 

1.58 

0.266004E+00 

1.59 

0,?570B4F+00 

1.60 

0 

.248437E+00 

1.61 

0.240058E+00 

1,62 

0.231946E+00 

1.63 

0.224095E+00 

1.64 

0.2165U1E+00 

1.65 

0 

•2091586*00 

1.66 

0.202061E+00 

1.67 

0.195204E+00 

1 .68 

0 . 18S5ROE+00 

1.69 

0.)82l84E+n0 

1.70 

0 

. 176009E + 00 

1.71 

0.17004BE+00 

1.72 

0. 164295E+00 

1.73 

0.15B743E+00 

1.74 

0.l53386E+nn 

1.75 

0 

.148218E*00 

1.76 

0.143231E+00 

1.77 

0.138420E+00 

1.78 

0.133779E+00 

1.79 

0.l29J02E*n0 

1.80 

0 

. 1249P2E*00 

1.81 

0.120815E+00 

1.82 

0.116794E+00 

1.83 

0.112915E+00 

1.84 

0.109171E+00 

1.85 

0 

.1055596+00 

1.66 

0.102073E+00 

1.87 

0.987078E-01 

1.88 

0.954596E-01 

1.89 

0.9232J9E-01 

1.90 

0 

.892963E-01 

1.91 

0.863727E-01 

1.92 

0.835671E-01 

1.93 

0.803973E-01 

1.94 

0.783554F-nl 

1.95 

n 

•759341E-01 

1.96 

0.736265E-01 

1.97 

0.714261E-01 

1.98 

0.693269E-01 

1.99 

0.673232E-01 

2.00 

0 

.654097E-01 

2.01 

0.635815E-01 

2.02 

0.618341E-01 

2,03 

0.601629E-01 

2.04 

0.585640E-01 

2.05 

0 

.5703346-01 

2.06 

0.555678E-01 

2.07 

0.541635E-01 

2.08 

0.528176E-01 

2.09 

0.515269E-01 

2.10 

0 

.5028P7E-01 

2.11 

0.491003E-01 

2.12 

0.479593E-01 

2.13 

0.468633E-01 

2.14 

0.4581UOE-nl 

2.15 

0 

.447975E-01 

2.16 

0.438237E-01 

2.17 

0.428B67E-01 

2.18 

0.419848E-01 

2.19 

0.411 164E-01 

2.20 

0 

.40279PE-01 

2.21 

0.394737E-01 

2.22 

0.386965E-01 

2.23 

0.379469E-01 

2.24 

0.372237E-01 

2.25 

0 

.365258E-01 

2.26 

0.358519E-01 

2.27 

0.352009E-01 

2.28 

0.345720E-01 

2.29 

0.339640E-01 

2.30 

0 

.333761E-01 

2.31 

0.328075E-01 

2.32 

0.322572E-01 

2.33 

0.317246E-01 

2.34 

0.312U8HE-01 

2.35 

0 

•307091E-01 

2.36 

0.302249E-01 

2.37 

0.297556E-01 

2.38 

0.293004E-01 

2.39 

0.288590E-01 

2.40 

0 

.284305E-01 

2.41 

0.2801476-01 

2.42 

0.276109E-01 

2.43 

0.272187E-01 

2.44 

0.26B375E-01 

2.85 

0 

.264671E-01 

2.46 

0.261069E-01 

2.47 

0.257566E-01 

2.48 

0.254158E-01 


TOTAL DRIFTS 0 . 557 T6f)F + 0 1 

TOTAL= 0.94A302E+U2 
BAKGNOl STOP 



RSP TEST CR3E 
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OSU / AARL 


AERIAL SPRAY PROGRAM 
LAST UPDATE 2/T/81 


ASP TEST case S 


program control 


N2D= 0 
NTUM= 0 

N30= 1 
NOIST= 1 

NEVAP= 1 
NPRINT= 0 

NPH0P2= 1 
NPLOT= 1 

NCM= 

EPS= 

1 

O.lOOOOOE-04 

parameter values 






A= 6.0000 

ZY0= 0.8000 

0= 0.0000 

CL= 0.6000 

ZZO= O.SOOO 

00= 1000.0000 

H= 6.0960 

DA= 0.11890IE + 01 

U= 53.0350 

VIS= 0.181330E-04 

g= 

VCW= 

9.8066 
1 .0000 

evaporation inputs 






pA= lonsa.io 

TOR= 23.B88900 

TWP= 18.333300 




PR0P2 INPUTS 






ZPR0P2= O.bOOU 

RP= 0.2000 

BK!= 1.7930 

PK?= 17.9300 

PK3= 0.0000 

PK4 = 

6.0290 

OISTRIRUTION INPUTS 






NCOL= 1 
NZNOZ= 1 
W0= 0.00000 

NROw= 6 
00=1000.00000 
SWInTH= 4.00000 

NDIAMN= 1 

x= o.oonoo 

UFRR= I. 00000 

NSTOEV= 1 
UD= 0.00000 

OwIOTH= 5.00000 

N0= 1 

VD= 0.00000 

nozzle no. tnoz 

ZNOZ 

DIAMN 

STOEV 0 


DRIFT 

1 -O.JOUO 

2 o.ouoo 

3 0.3000 

0.4000 

0.4000 

0.4000 

POO. 0(100 
200.0000 
PflO.OOOO 

50.0000 1.0000 
50.0000 1.0000 
50.0000 1.0000 


7.1442 

13.0321 

9.2220 



SINGLE droplet LATERAL DISPLACEMENT 


NO-I7LE! 

YG 

-0.320151E+00 
-0.'*07045E + 00 
-O.A59665e* 00 
-0.bb7H48E*00 
-0.7839A7E+00 
-U.IA702Ar+Ol 
-0. Ibl925f +01 
O.OOOOOOE-»Ol 
O.OOOOOOf+OI 
O.UOOOOOF+01 

o.ouooooe+uI 

O.UOOOOOe*U 1 
O.OOOOOOF+01 
0. OOOOOOE+Ol 
O.OOOOOOF+Ol 


1 

UIA 

0.6UOOnnE*03 
0. JUoOnOE+OT 
0.<i500onE*OT 
o.iUoOooE+ns 
0. lbOO(iOF + (n 
0. 1 uoOoOE + OT 
0,VM437SE+n? 
o.uMflOnnF+ni 
o.uuooooF+ni 
O.uuoonnE+o i 
'o.uoooooE+ni 

O.UUOOOOF+Ol 
0. ouoOonE + n i 
o.uui)OonE+ni 
o.uuoonoE+ni 


nlAG 

0.S98080F+03 

0.297088E+03 

0.246S94E+U3 

0.I9578UE+03 

0,I4A136E*03 

0.74o140E*o2 

O.G32314E+02 

0.000000t*0l 

O.OOOOUOE+Ol 

o.oooouoe+oi 

0.000008f* 0I 
O.OOOOUOE+Ol 
O.OOOOUUL+01 
0. OOOOOOE+Ol 

n.oooooOE+ui 


YQ 

-0.076441E-01 
-0.blS8b5E-0l 
-0,b39186e-01 
-0.593184E-01 
-0.842100E-01 
-0.139098E+00 
-0. 172918E+00 
0.J219O2E+00 
0.325157E+00 
0. 704499E-01 
0. OOOOOOE+Ol 
0. OOOOOOE+Ol 
0. OOOOOOE + Ol 
0. OOOOOOE + Ol 
0. OOOOOOE+Ol 


2 

DIA 

0.600000E+03 
o.3aooaOE+o3 
O.27bO0OE+O3 
0.2bOOOOE+03 
0.200000E+03 
0.15oOOOE+03 
0.12b0u0E+03 
O.IOOOOOE+03 
0.725000E+02 
0.O90625E+02 
0. OOOOOOE+Ol 
0. OOOOOOE+Ol 
0. OOOOOOE+Ol 
O.OOOOOOE+OI 
O.OOOOOOE+Ol 


DIAG 

0.b97780E+03 
0.296151E+03 
0.270796E+03 
0.245340E+03 
0.193927E+03 
0.140996E+03 
0.112875E+03 
0. n08b2E + 02 

0.461472E+02 
O.is7027E+O2 
0. OOOOOOE+Ol 
0. OOOOOOE+Ol 
0. OOOOOOE+Ol 
O.o'OOOOOE + 01 
0. OOOOOOE+Ol 


YG 

0.54474SE+00 
0.70fll99E + 0tl 
O.701975E+OO 
0.896006E+0D 
0.109198E+01 
0.152290E+01 
0t209856E+01 
0. OOOOOOE+Ol 
0. OOOOOOE+Ol 
0>OOOOOOE+01 
0. OOOOOOE+Ol 
0. OOOOOOE+Ol 
0. OOOOOOE+Ol 
0. OOOOOOE+Ol 
0. OOOOOOE+Ol 


3 

DIA 

0.600000E+03 
0.300000E+03 
0.250000E+03 
0.200000E+03 
0.150000E+03 
O.lOOOOOE+03 
0.770271E+02 
O.OOOOOOE+Ol 
0. OOOOOOE+Ol 
O.OOOOOOE+Ol 
O.OOOOOOE+Ol 
O.OOOOOOE+Ol 
O.OOOOOOE+Ol 
O.OOOOOOE+Ol 
O.OOOOOOE+01 


DIAG 

0.598063E+03 
0.297044E+03 
0.246526E+03 
0.196685E+03 
0.144039E+03 
0.900272E+02 
0.613805E+02 
O.OOOOOOe+ 01 
0. OOOOOOE+Ol 
O.OOOOOOE+Ol 
0. OOOOOOE+Ol 
O.OOOOOOE+Ol 
O.OOOOOOE+Ol 
O.OOOOOOE+Ol 
O.OOOOOOE+Ol 



final distribution 


YG 

conceniratton 

YG 

Concentration 

YG 

concentration 

YG 

concentration 

YG 

concentration 

-1.51 

0.919440E-02 

-1.50 

0. 103502E-01 

-1.49 

0.116264E-1)1 

-1.48 

0.1302B9E-01 

-1.47 

0.145629E-01 

-1.46 

0.162J40E-nl 

-1.45 

0.180510E-01 

-1.44 

0.200232E-01 

-1.43 

0.221601E-01 

-1.42 

0.244715E-01 

-1.41 

0.269678E-01 

-1.40 

0.296594E-01 

-1.39 

0.325572E-01 

-1.38 

0.356725E-01 

-1.37 

0.390157E-01 

-1. 16 

0.426016E-01 

-1.35 

0.464390E-01 

-1.34 

0.505412E-01 

-1.33 

0.549204E-01 

-1.32 

0.595890E-01 

-1.31 

0,645596E_nl 

-1.30 

0,698445E-01 

-1.29 

0.754561E-01 

-1.28 

0.814056E-01 

-1.27 

0.87707BE-01 

-1.26 

0.943714E-01 

-1.25 

U.10I408E400 

-1.24 

0.108B29E400 

-1.23 

0.116644E400 

-1.22 

0.124861E+00 

-1.21 

0.133489E+00 

-1.20 

0.142533E400 

-1.19 

0.152000E+00 

-1.18 

0.161892E+00 

-1.17 

0.172212E+00 

-1.16 

0.182959E*00 

-1,15 

0,1941318400 

-1.14 

0.205723E+00 

-1.13 

0,217726E+00 

-1.12 

0.230131E+00 

-1.11 

0.2429<;3E*00 

-1.10 

0.2560P5E400 

-1.09 

0.269595E+00 

-1.08 

0.2B342BE+00 

-1.07 

0.297556E+00 

-1.06 

0.31194bE+n0 

-1.05 

0.326556E400 

-1.04 

0.341349E+00 

-1.03 

0.356277E+00 

-1.02 

O.371209E+OO 

-1.01 

0.3863J1E+00 

-1.00 

0.401346E+00 

-0.99 

0.416273E+00 

-0.98 

0.431049E*00 

-0.97 

0.445610E+00 

-0.96 

fl,459892E+nO 

1 -0.95 

0.473831E400 

-0.94 

0.48736BE+00 

-0.93 

0.500447E+00 

-0.92 

0.513020E+00 

-0.91 

0.52504dE*nO 

-0.90 

U.536505E+00 

-0.89 

0.547382E+00 

-0.88 

n.557687E+00 

-0.87 

0.567455E+00 

-0.86 

0.576747E+00 

-0.85 

0.585659E+00 

-0.84 

0.594326E+00 

-0,83 

0.60292BE+00 

-0.82 

0.611698E+00 

-0.81 

0.6209J1E+00 

-0.80 

0.630991E400 

-0.79 

0.642328E+U0 

-0.78 

0.655514E+00 

-0.77 

0.671325E+00 

-0.76 

0.690124E+00 

-0.75 

0.712217E400 

-0.74 

0.737946E+00 

-0.73 

0.767698E+00 

-0.72 

0.801901E+00 

-0.71 

0.84103BE*00 

-0.70 

0.8a5644E*00 

-0.69 

0.936312E+00 

-0.68 

0.993696E+00 

-0.67 

0.105B52E+01 

-0.66 

0.113158E*01 

-0 . 65 

0.121374E401 

-0.64 

0.130592E+01 

-0.63 

0.140916E+01 

-0.6? 

0. 152453E+01 

-0.61 

0.165J20E+01 

-0.60 

0.179643E+01 

-0.59 

0.195549E+01 

-0.58 

0.213168E+01 

-0.57 

0.232615E+01 

-0.56 

0.253965E+01 

-0.55 

0.276672E+01 

-0.54 

0.299187E+01 

-0.53 

0.320669E+01 

-0.52 

0.340104E+01 

-0.51 

0.3562BlEt01 

-0.50 

0.367833F401 

-0.49 

0.373321E+01 

-0.48 

0.371375E+01 

-0.47 

0.360894E+01 

-0.46 

0.341<i94E + ol 

-0.45 

0.313619E401 

-0.44 

0.279241E+01 

-0.43 

0.238913E+01 

-0.42 

0.194197E+01 

-0.41 

0.147452E+01 

-0.40 

0.100968E401 

-0.39 

0.581835E+00 

-0.38 

0.260214E+00 

-0.37 

0.796776E-01 

-0.36 

0.136499E-01 

-0.35 

0.938853E-03 

-0.17 

0.168032E+01 

-0.16 

O.221909E+O1 

-0.15 

0.265141E+01 

-0.14 

0.302279E+01 

-0.13 

0.430308E401 

-0.12 

0.496400E+01 

-0.11 

0.561741E+01 

-0.10 

0.618396E+01 

-0.09 

0.657872E+01 

-0.08 

0. 146893E402 

-0.07 

0.128646E602 

-0.06 

0.969981E+01 

-0.05 

0.625907E-01 

-0.04 

0.601423E-nl 

-0.03 

0.577604E-01 

-0.02 

0.554479E-01 

-O.Ol 

0.532070E-01 

0.00 

0.510397E-01 

0.01 

0.4894^6E-nl 

0.02 

0.46931BE-01 

0.03 

0.449931E-01 

0.04 

0.431322E-01 

0.05 

0.413491E-01 

0.06 

0'.396439E_o1 

0.07 

0.380163E-01 

0.08 

0.365174E-01 

0.09 

0.350553E-01 

0.10 

0.336704E-01 

0.11 

0.323619E-nl 

0.12 

U.311288E-01 

0.13 

0.299698E-01 

0.14 

O.280839E-O1 

0,15 

0.278698E-01 

0.16 

0.269261F-nl 

0.17 

0.260517E-01 

0.18 

0.252451E-01 

0.19 

0.245051E-01 

0,20 

O.2303O7E-OI 

0.21 

0.232805E-IU 

0.22 

0.226735E-01 

0.23 

0.221889E-01 

0.24 

0.217659E-01 

0.25 

0.214037E-01 

0.26 

0.211019E-01 

0.27 

U.2U8601E-01 

0.28 

0.206784E-01 

0.29 

0.205569E-01 

0.30 

0.204959E-01 

0.31 

0.204961E-01 

0.32 

0.2U5586E-01 

0.61 

0.398603E-03 

0.62 

0.188731E-02 

0.63 

0.691297E-02 

0.64 

0.204101E-01 

0.65 

0.5027?8e-01 

0.66 

0.106330E+00 

0.67 

0.197834E+00 

0.68 

0.330412E+00 

0.69 

0.503837E+00 

0.70 

O.711515E400 

0.71 

0.941824E*00 

0.72 

0.11B517E+01 

0,73 

0.143373E+01 

0.74 

0. 167951E + 01 

0.75 

0.191594E401 

0.76 

0.213787E+01 

0,77 

0.234145E+01 

0.78 

O.252309E+O1 

0.79 

0.268062E+nl 

0.80 

0.280446E401 

0.81 

0.289454E401 

0.82 

0.295156E+01 

0.83 

0.297733E+01 

0.84 

0.297448E+nl 

0.85 

0.294617E+01 

0.86 

0.289584E+01 

0.87 

0.282701E+01 

0.88 

0.274313E+01 

0.89 

0.?64744E+nl 

0.90 

0.254306E401 

0.91 

0.243368E401 

0.92 

0.232165E+01 

0.93 

O.22O064E+O1 

0.94 

0.2096U5E+nl 

0.95 

0. 198500E401 

0.96 

0.18763BE+01 

0.97 

0.177087E+01 

0.98 

O.166097E+O1 

0.99 

0.1571U6E+01 

1.00 

0. 147738E + 01 

1 .01 

0.136806E+01 

1.02 

0.130318E+01 

1.03 

0.122273E+01 

1.04 

0.114665E*nl 

1.05 

0.1074B6E+01 

1.06 

0.100723E+01 

1.07 

0.943617E+00 

1.08 

0,eB3859E*00 

1.09 

0,827786t+n0 

1.10 

0.775062F400 

1.11 

0.725256E+00 

1.12 

0.678293E*00 

1.13 

0.634094E+00 

1.14 

0.8925b 5f4.no 

1.15 

0.55J596E*00 

1.16 

0.517082E+00 

1.17 

0.482905E+00 

1.18 

0.450947E+00 

1.19 

0.421089t4n0 

1.20 

0.393215E*00 

1.21 

0.367209E+00 

1.22 

0.342959E+00 

1.23 

0.320356E+00 

1.24 

0, 2992988400 

1.25 

0,?79683F400 

1 .26 

0.261418E+00 

1.27 

0.244412E+00 

1.28 

0.228580E+00 

1.29 

0.2138928400 

1.30 

0.200122E400 

1.31 

0.187350£*00 

1.32 

0.I75459E+00 

1.33 

0.16438BE+00 

1.34 

0.1540788400 

1.35 

0.144974E400 

1 .36 

0.135528E+00 

1.37 

0.127192E+00 

1,38 

0.119422E+00 

1.39 

0. 1 12n8f 400 

1 .40 

0.105422E400 

1.41 

0.991196E-U1 

1.42 

0.93P381E-01 

1.43 

0.877476E-01 

1 .44 

0.8262028-01 

1.45 

0. 77RJ01F-O1 

1.46 

0.733534E-01 

1.47 

0.691679E-01 

1.48 

0.652532E-01 

1,49 

0,6159028-01 

1 .50 

0.5H16I3F-01 

1.51 

0.549501E-01 

1.5? 

0.519417E-01 

1 .53 

0.4912S5E-01 

1 .54 

0,4649^38-01 

1.55 

U.490436E-01 

1.56 

0.417504E-U1 

1.57 

0.396055E-01 

1.58 

0.375977E-01 

1.59 

0.357167E-01 

1 .60 

0.3J9532E-O1 

1.61 

0.322984E-01 

1 .62 

0.307446E-01 

1.63 

0.292843E-01 

1 .64 

0.279110F-O1 

1 65 

0.266 1 B5F-01 

1 .66 

0.2S4011E-01 

1 .67 

n.?4?53BE-01 

1 .68 

0.231716E-01 




k • • * 


1 . < u 

•> . r 1 1 o- 1 

1 . • 1 

' 

1.74 

0.i7HaJ3E-nl 

1.7S 

<1.170913E-01 

1.76 


1.79 

0.l45au4E-fll 

1.80 

0.139563E-01 

1.81 


1.S4 

0.119615E-nl 

1.85 

0. llbao8E-01 

1.86 

( 

1.89 

o.995a.:5E-oa 

1.90 

0.960347E-0a 

1.91 


1.94 

0.8J5502E-0? 

i.-is 

O.8075B0E-Oa 

1.96 


1.99 

o.707lo4E-na 

?.00 

O.684512E-0a 

a.oi 

( 

a. 04 

0.6028a5E-0? 

a. 05 

0.S84368E-0a 

a. 06 


a. 09 

0.817J40E-O? 




( 

TOTAL DPIFT= 

0.97994JE+nl 





TOTALS 0.fl73970E*0<! 
BAKGND! STOP 


( 


( 

Ln 

^ ( 
C 
( 
( 
{ 
Iv 
( 
(. 


I 


b*<U£l*«UL.-Wi 

i • t c 


i • / J 

U* 


0.163973E-01 

1.77 

0.157389E-01 

1.78 

0.151139E-01 

■y 

0.134200E-01 

1.82 

0.129097E-01 

1.83 

0.124241E-01 


0.111006E-01 

1.87 

0.106998E-01 

1.88 

0.103174E-01 


0.927016E-oa 

1.92 

0.895150E-oa 

1.93 

0.864670E-02 

■i 

0.780837E-02 

1.97 

0.755215E-02 

1,98 

0.730654E-02 


0.662831E-02 

a. 02 

0.642017E-02 

a. 03 

0,622028E-oa 


0.566624E-0? 

2,07 

0.549559E-0a 

2.08 

0,53314lE-0a 

) 


) 

1 

) 

) 

) 

) 

) 

) 

) 

) 

) 

) 

) 

) 

) 

) 

) 
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APPENDIX B 


FORTRAN COMPUTER PROGRAM LISTING 
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o o o o 


c 10 

C ASP _ aerial spray program 20 

c 30 

C ASP WAS RFEN nESIGNEO TO PrpOICT THE GROUND DEPOSITION 40 

C «NO percent of material lost due to drift from an 50 

c Agricultural aircraft* the program is written to input ,so 

c umta in si units and includes the following options: 70 

C 1) ao OR 30 wake model and droplet dynamics 80 

C 2) a ORODLET EVAPORATION MODEL 90 

c 3 ) A propeller Slipstream model loo 

C 4) A CPOSSWINO MODEL 110 

C S) A tunnel wall model 120 

c G) ThF ARILITY FOR the USER TO SUPPLY HIS OWN 130 

C FLOWFIELD model. 140 

c 150 

C documentation available in the form of a NASA CONTRACTORS 160 

C RtPORT 170 

C 180 

C UOOE written at THE AERONAUTICAL AND AsTRONAUTICAL ENGINEERING 190 

C research LARORATORY* the OHIO STATE UNIVERSITY. COLUMBUS. 200 

C OHIO 19R1 author MICHAEL 9* BRAGG 210 

C 220 

C LAST update 2/1/81 230 

C 2A0 

C 250 

CDC PROMRAM mAIN(INPuT.0UTPUT»TAPE8.TAPE5=INPUT»TAPE6=0UTPUT) 260 

DIMtNSION YMOZ(IOO) . Sl Zt ( 1 5 ) .01 ST ( 1 0 0 1 ) , DR I FT ( 1 0 0 ) .NRUN U 0 0 ) 270 

DIMENSION DIAMN(IOO) .STOEVdOO) .0(100) .ZNOZ(IOO) .DIAGdS.lOO) 28 0 

DIMENSION IRUN(SO) .TRAJdS.lOO) .SiZdS.lOO) 290 

K .SLPdS) .XlIST (5) ,oliST(5) 300 

DIMENSION YY(13.13) .SAVE<13.13) .CSAVEd3.3) .YMAX(13) .MET (100) 310 

DIMENSION ERROR (13) . 0 Y d 3 ) , yY 1 (1 3 ) .PW ( 13 . 1 3 ) . I P d 3 ) .TITLE (13) 320 

COMMON /AREAI/ A,CL.8»U,G,2YO,ZZO.OA.VIS.VCW.O.DO 330 

COMMON /AREA2/ N20 .N3D » NE VAP . EPS . NOR0P2 »NCW . N TUN ♦ NO I ST . NPR I NT 340 

COMMON /4RFA3/ PS WR , P AE , P v . SCN . 0 V 350 

COMMON /areas/ 2PROP2.PK1 .PK2,PK3.PK4.RP 360 

COMMON /AREA6/ JJ . I RUN . I COUNT 370 

COMMON /AREA7/ NVAR 380 

COMMON /AREAIO/ TRAJ.DIAG.SIZ 390 

400 

INPUT THE NECESSARY PARAMETERS ACCORDING TO THE OPTIONS 410 

StLECTED and OUTPUT hEADER ANO TITLE 420 

430 

800 continue 440 

REWIND 8 450 

REAO(5.3.END=700) ( T I TLE ( K ) ,K^ 1 , 13 ) 460 

3 F0RMAT(13A6) 470 

wRIIt(6,S) (TITLE(K) *K=1»13) 480 

5 FORMAT(miii,t50."OSU / A ARL''/*'0” . T50 « 490 

h "AERIAL SPRAY PROGRAm"/"0"«T50» 500 

& "LAST UPDATE 2/1/81"/////"0"»13A6//) 510 

DO I Ksl.ion 520 

DRIET(K)=o. 530 

NRUN(K)aO 540 

DO I J=1.15 550 
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nan 


SIZ<vJ*K) = 0. 560 

0IAti(J.K)=0. 570 

1 TRAU» J.(f)=0. 580 

00 K=l,10Ol 590 

2 OISl (K)=0. 600 

ICoOnT=i 610 

JJ=U 620 

REAU <5f lOlNOPOPS 630 

10 FOPMATdS) ' 640 

00 410 I=1.nDR0PS 650 

REAO(5»?0) N2D.N30.NEVAP»NPR0P2.NC-<»NTUN,N0ISrt 660 

6 NPRINT*NPLOT»tPS 670 

20 F0RMAT(oip,?xtEl5.6) 680 

IFlNtVAD.EQ.l) READ(5*30) PA»TD0«TW0 690 

30 FORMAT nno. 6) 700 

IF (NPPOPP.EO.l ) PEAO(b*50) ZPPOP2tPKl ,PK2tPK3»PK4,PP 710 

50 FORMAT<6F10.4) 720 

REA0(5.60) A.CLtfl,U.G«ZYO,ZZO 730 

60 FORMAT(7F10.4) 740 

REAU(5.70) DAtVIS«VCW*0 750 

70 FORMAT(2E20.6*2F10.4) 760 

IF(NUISt.EO.I) go to 90 770 

REAU15»80> DIA,D0*X»Y*Z»UD»VD.W0 780 

80 F0RHAT{RFin,5) 790 

GO TO 120 800 

90 CONTINUE 810 

REAU(5«R5) MCOL*NROWfNOlAMN»NSTDEV»NQ«NZNOZ 820 

95 F0RMAT(6I5) 830 

REAIH5»ioO) Dn,XO,UDO*VDO*wDO*5WlDTH»OERRtOWlDTH 940 

100 FORMAT(flE10.5) 850 

REAU(5.110) (YN0Z<K) »K=1.NC0L) 860 

110 FORMAT(8F10.5) 870 

REAU15.U0) (SIZE (K) ,K=1 *NR0W) 880 

IF (NOIAmm.EO. 1 ) PEAO(S*110) (DIAMN (K) «K=1 fNCOL) 890 

IF(NUIAmn.Eq.O) REAO(b.llO) OIAMnI 900 

IF(NSTDfv.EO.I) PEAO(5»llO) (STOEV (K) «K=1 .NCCU 910 

IF(NSTOEV.EO.O) RFaD(5»110) STOEVI 920 

IF(NO.EO.I) REA(?(5*110) ( Q (K ) ,K = 1 .NCOL) 930 

IF(NO.EO.O) PEAO(5,110) 01 940 

IF(N4N07.EQ.1) READ(5»110) (ZNOZ(K) .K=1»NC0L) 950 

IF (NZNO7.FQ.0) READ(5»ll0) ZNOZl 960 

no U5 K=l.K'COL 970 

IF (NUIAmn.EO.O) DIAMN<K)=DIAMN1 980 

IF(nQ.E(3.0) 0 (K)=Q1 990 

IF(NsTOfv.EO.O) STDEV(K>=ST0EV1 1000 

IF (NZNOZ.EQ.O) ZNOZ(K)=ZNOZl 1010 

115 CONTINIJf 1020 

120 continue 1030 

C=3»i4l5Q26935a9793 1040 

1050 

Calculate evaporation parameters 1060 

1070 

IF (NtVAP.EQ.O) 60 TO 130 1080 

T0Ht< = TDa*27->, 15 1090 

T0BM=TDR*(9./5. ) *3?.+A59.67 1100 
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TWBH=TWR*(9./5.)+32.+459.67 UlO 

PA£=PA*1 .4504E-4 1120 

PS=tAP(54.6-?2O-ia301 .6atl/TD9R-5.16923*AL0G(T0BH) ) 1130 

PS-(O=exp(54,f,329-123 0l.688/TWnR-5.16923*AL0G(TrtBR) ) 1140 

PVaHbw8-( t.?405*(PAE-PSW9) )/(. 621 94* ( 1 0 75 . 8963- . 56983* 1150 

<TWBR-491,69) ) ) )*(TUBR-TWBR) 1160 

RH=PV/P<; 1170 

DA=( (PAE*144,)/(g*3. 260839*53. 34*T0“R) )*515.4 1180 

DV=«b.2PE-G*TDBK**l .88) *1 ,E-4 1190 

vIS=1.795E-5*(T0RK/293.16) **1.5* ( ( 293 . 16 + 1 1 0 .)/( T08K + 1 1 0 . ) ) 1200 

SCN=10V«0A) /VIS 1210 

130 continue 1220 

1230 

OtTERMINE NVAR. ThE NUMBER OF SIMULTANEOUS FIRST ORDER 1240 

differential equations to re SOLVED 1250 

1260 

MVAR=4 1270 

IF(N30.nE.O) NVAR=6 1280 

IF(NtVAp,EO. 1 ) NVAR=NVAR+1 1290 

IF«N2D.fQ. 2.0R.N3D.EQ.2) GO To 138 1300 

IF (NTUN.eQ . 1 ) GO TO 136 1310 

IF(NCW.eo.1.AND.NPROP2.EQ.1) go TO 132 1320 

IF(NLW.EQ.I) GO TO 134 1330 

IF(NPROPP.EQ.I) NVAR=NVAR+1 1340 

GO TO 138 1350 

132 NVAR=NVAR+6 1360 

GO TO 138 1370 

134 NVAR=NVAR+4 1380 

GO TO 138 1390 

136 NVAH=NVAR+2 1400 

IF<NRR0P2.EQ.1) NVAR=NV«R+1 1410 

GO TO 13P 1420 

138 CONTINUE 1430 

1440 

OUTPijT THE constants THAT WERE INPUT 1450 

1460 

WRpt(6,l50) N2D.N3D.NEVAP.NPR0P2tNCW*NTUN.N0IST.NPRINT« 1470 

t. NPL0T.EPS 1480 

ISO FORMAT(>io''. "PROGRAM CONTROL".//" " . T5, "N2D=" . 1 2 . T25 , "N3D=". 12 » 1490 

i t45."NEVAP=".I2.T65."Npr0P2="»I2.T85."NCw=".I2/" " 1500 

& .T5."NTUN=".I2.T25."NDIst=".I2»T45."NPRINT=". I2.T65. 1510 

t "NPL0T=".I2.T85."epS=",E13.6/> 1520 

wRIiE(6.l60) A.CL»B*U*6*ZYO.ZZO.DA.ViS.VCW.O.OD 1530 

160 FORMAT("0"."P4RaMETER VALUES"//" "«T5,"A=",F10.4.T25."CL="» 1540 

A F10.4,T45."B="»F10.4.T65."U="»F10.4.T85."G=".F10.4/" 1550 

S. T5."ZY0=".F10.4.T25."ZZn=".Fl0.4,T45."OA=",E13.6.T6b, 1560 

A "VIS=i>,E13.6.T85,"VCW=",E10.4/" " , T5 . "U=" ,F 1 0 . 4 , 1570 

A T25."DD=",F10.4/) 1580 

iF(NtVAp.EQ.l) WRITE(6 i170) PA.TDS.TWB 1590 

170 format ("0". "EVAPORATION InPUTs"//" " » T5» "P A=" .FI 0 . 2. T25. 1600 

A "tD8=",F10.6.tA5."TW9s:",f10 ,6/) 1610 

IF(nPR0p2.E0.1) wRITE«6»190) 7PROP2.PK 1 .PK2 .PK3»Pk4 , rp 1620 

190 format ("n"."PROP2 INPUTS"//" ", T5» "ZPROP2=" » Fl 0 .4, T25. 1630 

A "DK1=",F10.4,T45.»PK2=",F10.4.T65,"PK3=",F10.4 1640 

A ,T85."PK4=".F10.4/II "»T5."RP=".F10,4/) 1650 
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c UIANOE FPOM rectangular TO ELLIPTICAL LOAD DISTRIBUTION 1660 

C UT AOJUSTING CL 1670 

CL= (‘♦./O *CL 1680 

IF(NDIST.EO.I) go to 320 1690 

wRITt(6,19S) X.Y,Z«UOtVOiWD.DlA 1700 

195 format (iiOH.itINITlAL vALUtS''//ii «,T5 . ''X=",F1 0 . 6 * T35 « 1710 

i ••y=".F10,6»TA5«'<Z="«F1o.6.T65*'*UD='**F10.6*T85» 1720 

i "VD=",F10.6/" "«t5«"W0='',F10.6.T25."DIA=», 1730 

^ F10.4/) 1740 

C 1750 

C CALL 9UB2D OR SUB30 TO CALCULATE THE DROPLET TRAJECTORY 1760 

C 1770 

IF(N2D.NE.0) call SUB2D<NVAR,Y,Z*V0.WD,T,DIA,YY,SAVE 1780 

& ,CSAVE,YMAx*ERROR.OY,YY1»PW«IP) 1790 

IF(N30.nE.O) call SU830(NVAR,X.Y.Z,U0,VD,WD.T*0IAtYY,SAVE 1800 

.CSAVE.YMAA, ERROR, 0Y«YY1*PW. IP) 1810 

lC0UNT=rC0UNT*l 1820 

WRITE(6,200) X* Y.Z»UD«VO»wO»Oia*T 1830 

200 format ( itQ". "FINAL VALUEs"//" " ,T5« "X=" * FI 0 . 6 ♦ T35 * " Y=" * 1840 

«. Fl0.6*T45.''Z="«Fl0.6,T65."UD="»F10.6»T85t"VD='** 1850 

t, F10.6/" ",T5."WO=".F10.6«T25*«DIA=",F10.4,T45. 1860 

f. "T=",F10,6/) 1870 

210 CONTINUE 1880 

GO TO 320 1890 

220 CONTINIJF 1900 

C 1910 

C this part OF THE MAIN PROGRAM CONTROLS THE NOISTsl CASE 1920 

C 1930 

WRI rt(6,22«5) NCOL*NROW»NDIAMN,NSTOEV*N(J.NZNOZiOD» 1940 

(, ' XOiUOO< VOO,wOO»SWIDTH.OERR*OWIDTH 1950 

225 format (" n". "DISTRIBUTION INPUTS”//” ”*T5»”NC0L=”. 12* I960 

«, T25.”MR0Wa”*I2*T^5,"NDlAMN=”*I2*T65*”NSTDEV=",I2. 1970 

% TP5."NQ=",I2/” "*T5*"N2NOZ=”«I2*T25*"DD=”*F10.5« 1980 

i T45*"Xa”,Fl0.5*T65t”UD="tFl0.5,T85*”VD=”*F10.5/ 1990 

i ” ”*T5."wn="*Fl0.5*T25,"SWlOTH="*F10.5«T45."DERR=”* 2000 

i F10.5.T65*”OWIOTH=",F10.5) 2010 

C 2020 

C HtRFoRM all the TRAJECTORIES REQUIRED TO CALCULATE THE 2030 

C HEQUIPEO DISTRIBUTION 2040 

C 2050 

DO CSO I=l,MCOL 2060 

DO 240 j=1,nROW 2070 

0IA=SIZE(J) 2080 

OIA1=OIa 2090 

K=K-1 2100 

T=0. 2110 

X=XU 2120 

UD=UD0 2130 

VD=YD0 2140 

WD="D0 2150 

Z=ZNUZ(n 2160 

Y=YNUZ(I) 2170 

IF<N20.NF.O) CALI. SU02D(NVAR,Y,Z*VO.«<D.T.OIA,yY,SAVEf 2180 

f. CSAVE.YMAX, ERROR, DY«YY1»PW, IP) 2190 

IF(N3D.NF,0) CALL SUB30(NVAH*X.Y*Z.U0.VD*WD*T*0IA*YY*SAVE* 2200 
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& CSAVE.YMAX»ERR0R,DY.YY1.PW*IP) 2210 

TRAJ(J«I)=y 2220 

Ol A<5 4 J» I ) =DT A 2230 

SI2l«J*I)=SIZE(J) 2240 

IF(Z.NE.o.) (50 TO 241 2250 

NRUN( I)=^RU^J( I) +1 2260 

IF (Z.EO.O.AnD.ABS (Y) .GT.DwIOTH) GO TO 241 2270 

240 CONTlNUF 2280 

C 2290 

C THIS SECTION UP TO STMT NO 250 CONTROLS THE CALCULATION 2300 

C Oh additional trajectories in order to estimate the drift 2310 

C «ND inserts particles in regions of a large gradient 2320 

c UUIA/dYG 2330 

c 2340 

KMISS=0 2350 

SIZMIS=n. 2360 

241 CONTINUE 2370 

J=NRUN(I)+1 2380 

METII)=0 2390 

IF(J.EO.l) nPIFT(I)=SlZ(J,I) 2400 

IF(J.EQ.l) Of) TO 250 2410 

OSIZl = Aqs(DlAMN(I)-SlZ(J-l*D) 2420 

K250=0 2430 

GO TO 243 2440 

234 CONTINUF 2450 

IF(Z.EO.o.ANO.ARS (Y) .OT.OWIDTH) GO TO 250 2460 

IF( (0SI71/STDEV(I) ) .LT.4M GO TO 242 2470 

ORIFTdlaO. 2480 

GO TO 2sn 2490 

242 CONfiNUE 2500 

C 2510 

C define acceptable ERROR IN TERMS OF STDEV 2520 

C 2530 

ERH=*S*STOEV(I) 2540 

IF(OSIZ1.LT.(2.*STOEV(I))) ERR=.30»STDEV < I ) 2550 

IF(USIZ1.LT.STDEV(I) ) EHR=.20ostDEV ( I ) 2560 

ERR=tRR*OERP 2570 

246 CONTINUE 2580 

IFCZtNE.n.) KMISS=1 2590 

IFIZ.NE.O.) SIZMIS=SIZ(J*I) 2600 

IF(J.GE,3) go to 239 2610 

IF(J.EQ.2) SIZ(J.I)=SlZ(J-l,I)-(SlZ(J-l.I)-SIZMIS)/2. 2620 

GO TO 245 2630 

239 CONTINUE 2640 

C 2650 

C UETEpmINF particle size for DRIFT ESTIMATE 2660 

C and calculate TRAJECTORIES 2670 

C 2680 

YGM1=TRAJ ( J-2, I ) 2690 

yG=TRAJ(j-1,i) 2700 

Yl=AbS(l./(YG-YNOZ(I) ) 1 2710 

Y2=ABS<1 ./(yGMI-YNOZ(I) ) ) 2720 

SLOh't3<siZ(j_2,n-SIZ<J-ltI) )/(Y2-Yl) 2730 

SlZtJ,I)=si7(J-l,I)-SL0PE*Yl/2 2740 

IF( tSL0pP*Yl/2.) .LT.ERR) SIZ(J.I)=SIZ(J-1.I)-ERR 2750 
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IF(blZ(j.I) .GT.SIZHIS) GO TO 245 2760 

SIZld«I)=SIZMiS*(SIZ(J-l.I)-SlZMIS)/2. 2770 

245 CONtiM'JE 2700 

DIA=SIZ(J,I) 2790 

OIA1=OIa 2000 

T=0. 2810 

X=XO 2820 

UD=UO0 2830 

VD=VD0 2040 

WD=*«D0 2850 

Z=ZNUZ(I) 2860 

Y=YN0Z(T) 2870 

IF (N20.NF-.0) CALL SU82D<NVAR,Y,Z*VD*rtO,T,DIA,YY,SAVe, 2830 

6 CSAvE,YMAA, ERROR, DY tYYl«PW*IP) 2890 

IF(N30.NE.O) call SUB30(NVAR»x.Y,Z»UD,VD»WD,T*0IA»YY,SAVE, 2900 

(. CS4VE,YMAX, error, 0Y,YY1,PW»IP) 2910 

2920 

STORf trajectory information, if ERR0R.lt. err go to 250 2930 

OTHERWISE ITERATE 2940 

2950 

TRA0(J,I)=Y 2960 

DIAO<J,I)=OIA 2970 

IFCMtTln ,EQ.l) GO TO 228 2980 

IF{4.NE.n.) GO TO 246 2990 

^RUN^I)sMRUN{I)+l 3000 

IF(ABS< Y) .LT.OWIDTH) GO To 330 3010 

GO TO 200 3020 

330 continue 3030 

IF( «SIZ(J,I)-SIZMIS) .GT.ERR) go TO 244 3040 

DRI^T<I)=SIZ(J,n 3050 

K25U=1 3060 

GO TO 243 3070 

344 continue 3080 

J=J+1 3090 

IF(J.LE.15) go to 246 3100 

WRITti6,?47) 3110 

247 F0RMAT(« ii.it*** max NUMBER OF RUNS EXCEEDED BEFORE DRIFT", 3120 

i " calculation COMPLETED") 3130 

ORI)-T(I)=SIZ(J,I) 3140 

60 fu 200 3150 

243 CONTINUE 3160 

3170 

THIS SECTION HANOLLES THE CASE WHERE OYG/DDIA CHANGES 3180 

SIGN RETWFEN points 3190 

3200 

IF(J.EQ.2) 60 TO 227 3210 

j2=u-2 3220 

no 410 K=l.j2 3230 

SLPlK> = (TRAj(K + l,I)-TRAJ(K,I))/(SlZ(K*l,n-SlZ(K,in 3240 

410 CONTINUE 3250 

DO 221 K=2.j2 3260 

IF(SLP<K) .GT.O.AND.SLPfX-l) .LT.O.) GO TO 222 3270 

IF(SLP(K) .LT.O.ANO.SlP tX-1) .GT.O.) GO TO 222 3200 

221 CONTlNUr 3290 

IF(Kc50,eQ.1> go to zSO 3300 
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GO TO 234 3310 

222 continue: 3320 

3330 

insert particles where DyG/dOIA changes sign 3340 

3350 

mETU) = i 33ft0 

IF(K250.PO. 1 ) GO TO 250 3370 

IADD=0 3380 

00 <?=2,J2 3390 

3400 

IF(SLP(k, .ut.O.AND.SLP(K-I) .LT.O.) GO TO 223 3410 

IF(iLP(K) .GT.O.AND.SLP(K-I) .GT.O.) GO TO 223 3420 

IFCIADD.nE.O) k=i<2 + IA00 3430 

DeLTA=Ans(SI7(K,I)-DlAMN(I) )/? tDEV(I) 3440 

ER=1.0«stDEV(I) 3450 

IF(0tLTa.LT.2.) ER=.75*STDEV(I) 3460 

IF(UELTa.LT.i .) ER=.50<‘STDEV<I) 3470 

ER=tH*DERR 3480 

IDELTA=4RS(SIZ(K+1.I)-SIZ(K,I) )/eR 3490 

IADD=IA004inELTA 3500 

OELTA=(8iz(K+l,i)-slZ<KtI) )/( idELTA+1) 3510 

Kl=K+l 3520 

Jl=>J-l 3530 

DO CiW k5=K1,J1 3540 

KK=K1-K5+J1 3550 

TRAJC<K4.I0ELTA, I)=TRAJ(KK,I ) 3560 

dial (KK*iOELTA* I )=DIAG(KK,I ) 3570 

SiZti'K+IOELTA, n=SIZ(KKtI) 3500 

224 CONTINUE 3590 

DO kk=1,IDELTA 3600 

SIZlK+KK,n=SIZ(KtI)+«K*OELTA 3610 

DIA=SIZ(K+KK,I) 3620 

T=0« 3630 

3640 

UO=UO0 3650 

VO=vD0 3660 

WD=wOO 3670 

2=ZNOZ(l) 3680 

y=YNOZ(I) 3690 

IF(N«;D.nE.O) call SUB^D (NVAR,y,Z.VO»WO,T*OIA.YYiSAVE. 3700 

K CSAVE.YMAX*ERROR,DY*YYlfPW.IO) 3710 

IF<NJD.np,0) call SUB30(NVAR,x,Y«Z*U0»VD»WD»T»DIA*YY*SAVE* 3720’ 

^ CSAVE*YMAx*ERROR,OYtYYl ,PW, IP) 3730 

TRAJ<KK+i<,I)=y 3740 

DIAG iKK + x, I ) =OIA 3750 

NRUN( I) =mrun< I) +1 3760 

3770 

226 continue 3780 

NRR=NRun(I) 3790 

223 continue 3800 

227 continue 3810 

3820 

ULTErmINE particle size for drift estimate and 3830 

CALCULATE THE TRAJECTORY 3840 

3850 
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DSIZi=Aq?(DTAMN(I)-SlZ(J-l#I) ) 3860 

IF ( (USI71/ST0EV ( I) ) .LT.4) go TO 237 3870 

ORII-TdisO. 3880 

GO TO 2Gn 3890 

337 continue 3900 

ERR=«5*STDEV ( I ) 3910 

IF(DSIZ1.LT. (2.»STDEV<I) 1 ) ERR=.30*STOEV(I) 3930 

IF(ObIZl.LT.STOEV(I) ) ERR=.20*STOEV ( I ) 3930 

ERR=tRR*nERR 3940 

336 CONIINUF 3950 

IF(^MISs.EO.n) 9I2(J.I)=SIZ(J-1.I)-1.1*ERR 3960 

IF(KMlS9.E0.n SIZ(J.I)=SIZ(J-l,I)-(SIZ(J-l*I)-SIZMIS)/2. 3970 

GO TO 348 3980 

328 CONTINUF 3990 

4000 

It: Particle within err go to 25o otherwise iterate 4010 

4030 

IF(Z.£0.0.) GO TO 239 4030 

KMIiS=l 4040 

SIZMIS=<;IZ(J. I) 4050 

GO TO 337 4060 

329 CONTINUE 4070 

NRUN(I)=nRUN(I)+1 4080 

IF(ARs(y) .LT.OWIOTH) eo To 340 4090 

ORIF f (Il=SIZ(J-l «I)-Aas( <SIZ(j,n-SlZtJ-l*I) )/(Y-TRAJ(J-It n > )* 410 0 

t. (0wIDTH-A8S(TRAJ(J-1,i) ) ) 4110 

GO TO 380 4130 

340 CONTINUE 4130 

IF(A8S(8IZ(J.n-SlZ(J-l*I)).GT,ERR) GO TO 233 4140 

0RIFT(I)=SI7(J*I) 4150 

GO To 280 4160 

333 CONTINUE 4170 

J=J*I 4180 

IF(J.LE.15) GO TO 236 4190 

wRITE(6,347) 4200 

250 continue 4210 

4220 

combine trajectory information to build-up the distribution 4230 

ANO estimate the drift output results 4240 

4250 

CALL swath (NCOL«NROW,UiOIAMN,STOEV,YNOZ»DWIDTH»NRUN»MET* 4260 

& SWIDTHtDISTiORIFT) 4270 

WRIIt(6,252) 4280 

252 FORMAT (MOM, iinoZZLE NO* "» 8x t "YN0Z”» 1 "ZNOZ" *8X * 4290 

i MniAMNM,i0X,"ST0Ev"*13x,"O''*12X, "DRIFT"/) 4300 

DO <i55 K=l«NCOL 4310 

wRITt (6,353) k.YnOZ(K) »ZN 0Z(K) ,0IAMN(K) tSTDEV (K) *Q (K) t DRIFT (K) 4320 

253 FORMAT(m ",4x,n*6F15.4) 4330 

?55 CONUNUF 4340 

WRITE(6,256) 4350 

256 FORMAT(mi"//" single DROPLET LATERAL DISPLACEMENT"//) 4360 

4370 

OUTPUT trajectory information to PRINTER 4380 

4390 

NP=NC0L/3 4400 
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NLEI" T=NcnL-NP*3 4410 

KP=-c 4420 

IFINP.EO.O) r,0 TO 425 4430 

OO 4U5 K1 = 1,NP 4440 ■» 

KK=1FIX(K1/3.)*2 4450 

IF(Kl.EO.KK) WRITE(6.415) 4460 

415 format (M 1 II * II II) 4470 

KP=Ki*3-2 4480 

KP1=AP+i 4490 

KP2=^°+2 4500 

lF(NtV4p.EQ.l ) 00 TO 450 4510 

WRlTt(6.420) KP,KP1.KP2 4520 

420 FORMAT{iioii,iiWOZZLE;»»Tl3,I3,T55»l3»T97,I3/" " • T7 , » Y6" * T2 1 ♦ "D I A" , 4530 

S. T49.iiYGi'.T63»"0IA'i,T91 ,iiYG"«TI05.i'OIA"> 4540 

00 4J0 K2=1.15 4550 

WRI It (6,440) TRAJ(K2.KP) ♦SIZIKE.KP) *TRAJ(K2«KP+1> . S I Z < K2 • KP+ 1 ) • 4560 

& TRAJ(K2.KP+2) ,SIZ(K2»KP+2) 4570 

440 FORMAT<ii "♦2F14,6.14X*2E14.6*14X.2E14.6) 4580 

430 CONTINUF 4590 

GO TO 405 4600 

450 continue 4610 

WRITt(6,460> KP.KP1.KP2 4620 

460 FORMAT(iio".i'NOZZLE!".T20»I3.T62.I3»T104.I3/'' "» T7 . "yG" * T2 1 » "DI A«, 4630 

S. T35.'iDIAGii.T49»"IG'i.T63."OIA".T77."DIAG".T91,i'YG"»T105. 464 0 

i "DTAh.T119.''DIAG") 4650 

DO “70 K2=1.15 4660 

WRITt(6,480) TRAJ(K2.KP) ♦ SI Z (K2 tKP) *0 1 AG (K2 .KP) »TRA J (K2 »KP*1 ) ♦ 4670 

i SIZ(K2*KP+1) .DIAG(K2»KP*1) *TRAJ(K2.KP42) »SIZ(K2»KP+2) 4680 

i .0IAG(K2*KP+2) 4690 

480 FORMAT(ii 'i,0E14.6) 4700 

470 CONIINUE 4710 

405 continue 4720 

425 CONTINUE 4730 

IF(NLEFt.EQ.O) go to bOO 4740 

KK=1FIX{K1/2.)*2 4750 

KP3=i'P + 3 4760 

KP4=KP*4 4770 

IF(KI.NE.KK) WRITE(6.415) 4780 

IFCNUEFt.EO.I ) GO TO 550 4790 

IF (NtVAP.EQ, 1 ) GO TO 510 4900 

WRIIE(6,490) KP3.KP4 4810 

490 format ( noil. iiMOZZLE! '*.T13. I3.T55. 13/" "»T7,"YG"« T21 . "01 A". T49»"YG" » 4620 

6 T63."DIA") 4830 

DO POO K?=l,i5 4840 

WRI IE (6. 440) TRAJ(K2.KP*3) , S I Z ( K2* KP+3 ) , TR A J ( K2 . KP+4 ) » SI Z (K2 .KP44 > 4850 

500 continue 4860 

GO To 6nn 4870 

510 continue 4880 

WRITt(6.520) KP3.KP4 4890 

520 FORMAT(iio", "NOZZLE)". T20*I3.T62.I3/ii " . T7»" Y6" »T21 . "DI A" . T35 » 4900 

6 "niAG".T49."YG"«T63."DlA".T77."DIAG") 4910 

00 530 K2=l,i5 4920 

WRI II t(6,480) TRAJ(K2.KP+3) ,SIZ(K2*KP43).DIAG(K2.KP+3)» 4930 

X TRAJ(K2,KP+4) ,SIZ(K2.KP+4) .0IAG(K2.KP+4) 4940 

530 continue 4950 
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GO TO 6no 4960 

550 continue 4970 

IF (NEVAp.eQ. 1 ) GO TO 580 4980 

WRITt(6,560) KP3 4990 

560 FOR«AT(mo».hmozzlE:h,T13»I3/ii n, t 7 » "YG" t T2 1 ♦ "0 1 A") 5000 

WHlTt(6,570) (TRAJ(K2fKp+3 ) ,siz(K2»KP+3) «K2=1*15) 5010 

570 FORHATIn ''t2El4,6) 5020 

GO to 600 5030 

580 continue 5040 

KRITt(6,585) KP3 5050 

585 FORM8T(mo''.»NOZZLE:"»T20»I3/" " , t7 »"YG" » T2 1 1 "D I A" . T35 1 "0 1 AG" ) 5060 

ViRITt (6,590) (TRAJ(K2*KP+3 ) .S I z ( K 2 f KP+3 ) ,DIAG (K2,KP+3) .K2=l i 15) 5070 

590 FORM4T(» ".3F14.6) 5080 

600 continue 5090 

WRIIE(6,260) 5100 

260 F0RM«T(him,»fINAL DISTRI0UTION"//"O"t5(3X*" Y6",6X. 5110 

6 "CONCENTRATION")/') 5120 

5130 

WRITE TO FILF 8 and OUTPUT TO 6 THE FINAL DISTRIBUTION 5140 

5150 

KJ=U 5160 

TOT=0. 5170 

START=-SWI0TH 5180 

NPTb=SWioTH*200*l 5190 

JJ=NPTS 5200 

DO 910 J=1,nPTS 5210 

XPU=STArt+(J-1)*.01 5220 

WRIlt(8,302) XPL.OIST(J) 5230 

302 F0RMAT(2E13.6) 5240 

TOT=tOT4.nIST< J) 5250 

lF(UlST(j) . lT.I.E-4) go to 310 5260 

IF(KJ.E0.5) KJ=0 5270 

KJ=KJ+1 5280 

XLIST (KJ)=XPL 5290 

OLIST(Kj)=OIST(J) 5300 

IFtKO.E-j.S) wrITE( 6.300) (XLIST (K) f DLIST (K) .K=1 t5) 5310 

300 FORMAT(h ",5(3X.F6.2t3x»E13.6) ) 5320 

310 continue 5330 

lF(Kd.NE.5) WRITEf6«300) (XLIST (K) .DLIST (K) «K=1 *KJ) 5340 

5350 

utTERMlNE TOTAL DRIFT 5360 

5370 

QF=U. 5380 

ORIFTF=o. 5390 

DO CPI K=l.NCOL 5400 

DRIt-TF=nPIFTF+ORIFT(K)«Q(K) 5410 

0F=QF+O(k) 5420 

261 continue 5430 

DRIhTF=OPIFTF/OF 5440 

WRI1£(6,?65> DRIFTF 5450 

TOT=TOT/oF 5460 

2g5 FORMAT(«0". "TOTAL DRIFT=".E13,6) 5470 

WRITL(6,266) tot 5480 

266 FORMAT (MO", ••T0TAL = "»EI^.P) 5490 

320 continue 5500 
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IF(NPlOt.EQ.I) call PUOTl (TITLE) 5510 

(30 TO 800 5520 

700 CONtiNUF 5530 

STOP 5540 

ENO 5550 

5560 

bURRoijTINE SUB2DS CALLED 3 y THE MAIN PROGRAM, SU820 CONTROLS 5570 

the trajectory calculation FOR ALL 20 CASES. GIVEN THE 5530 

Initial values and user selected options, sueao calculates 5590 

the trajectory and returns the final values of the variables 5600 

5610 

subroutine SUB20(NVAR,Y,Z,VD,wD,T,0IA,YY,SAVE,CSAVE, 5820 

«, YMAX*ERR0R,0y,YY1,PW,IP) 5630 

dimension YY(13,NVAR) »SaVE(13,NVAR) ,CSAVE (NVAR.3) ,YMAX(NVAR) 56A0 

DIMENSION ERDOR(NVAR) *OY(NVAR) ,YYl (NVAR) 5650 

DIMENSION PW(NVAR,NVAR) ♦IP(NVAR) ,IRUN(50) 5660 

COMMON /AREAl/ A.CL,BtU*G,ZYO,ZZO,OA,VIS,VCW,D,OD 5670 

COMMON /AREA?/ NED . N30 ,NE V AP , FPS ,NPROP2 , NCW , NTUN ,ND I ST ,NPRINT 5680 

COMMON /APEA5/ ZOROP2*PK1 ,PK2,PK3*PKA,RP 5690 

COMMON /AREA6/ JJ, IRUN, ICoUNT 5700 

COMMON /AREA8/ OIAl 5710 

5720 

initialize the dependent variables and STORE IN array 5730 

YY(I.nVAR) 57A0 

5750 

OIA=DIA»i .E-6 5760 

OIA1=OIa 5770 

C=3* IA15Q2653509793 5780 

T=0.E0 ,, 5790 

YY(1,1)=y 5800 

YY(i»2)=vD 5810 

YY(1,3)=7 5820 

YY(1,A)=wD 5830 

IF<NEVAp.E0.1 ) YY(1,NVAR)=0IA 5840 

IF(N2D.FQ.2» go to 4 5850 

IF(NVAR.PE.IO) GO TO 1 5860 

IF(NVAR,GE.8.AND.NPROPE«EQ.O) go to 2 5870 

IF(NTUN.EQ.1.AND.NPR0P2.EQ.1) YY(1,7)=ZPR0P2 5880 

iF(NfUN.FQ.l) GO TO 3 5890 

IF(NHR0P2.E0.1) YY(1,5)=ZpR0PP 5900 

GO 10 4 5910 

1 YY(I* 10)=ZPPOP2 5920 

YY(i*9)=0. 5930 

2 YY(1«8)=7Z0 5940 

YY(1,7)=-ZY0 5950 

3 YY(1*6)=7Z0 5960 

YY(1,5)=zY0 5970 

4 CONTINUf 5980 

1F(VU.EO.O.AND.WO.EQ.O.) call INC0N(YY,0IA,NVAR) 5990 

lF(YYtl,4 ) ,EO.O.EO.AND,yY(1,2) .EQ.O.EOJGO TO 20 6000 

6010 

initialize values EOR DIFSUB 6020 

6030 

H=l*E-4 6040 

MF=2 6050 
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00 b I = X.NVAR 6 O 1 SO 

YMAX(I)=1.E0 60T0 

5 continue 6080 

MAX0tR=#i 6090 

HMIN=1.e-16 6100 

HMAX=1.P0 6110 

JSTART=0 6120 

L00^'=0 6130 

YPROK=0. 6140 

6150 

HtRFoPM the step integration BY SUCCESSIVE CALLS OF DIFSUB 6160 

6170 

OIFSUR IS A FORTUOI SUBROUTINE TO SOLVE A SYSTEM OF N SIMULTaNE 6180 

OIFFERENTIAL equations 6190 

6200 

DO I'* K=],moO 6210 

DO 9 I=],NVAR 6220 

YY1(I)=yY(1,I) 6230 

9 CONIINUE 6240 

Tl=f 6250 

6260 

CALL 0IFSBM(MVAR,T, YYtSAVE*CSAVE»H»HMIN»HMAX.EPS»MF,YMAX.ERRORt 6270 

^ KFLAO* JSTART»MAXOep,PW» IP) 6280 

6290 

monitor VARIABLES AFTER EACH TIME-STEP TO CONTROL OUTPUT 6300 

AND STORAGE OF TRAJEOTORY AnO TERMINATE TRAJECTORY IF 6310 

particle PECOmES entrained in the VORTEX OR HITS THE GROUND 6320 

6330 

DIAPL=YY(l.MVAR)/l.E-b 6340 

iF(NfcVAD.EQ.l) YY(1.NVAR)=YY(1,NVAR)/1.E-6 6350 

IF(NHRINT.EO.I.AND.NDIST.EQ.O) WRITE (6.11) K . T. ( Y Y ( 1 . I ) . I=l .NV aR) 6360 

11 FORHAT(ii ••,I4.4X.12F10.6) 6370 

iF(NtVAp.EQ.l) YY(1 »NVAr)=YY(1,NVAR)*1.E-6 6380 

L=K-^JJ 6390 

IF(YY{1,3).lt.O.EO)GO TO 16 6400 

IF<NUIST.E0.1) go to 10 6410 

WRirt(8,22) T.YY(l.l) »YY(1,3) 6420 

22 F0RMAT(3F13,6) 6430 

10 continue 6440 

lF(NtVAp.E0.1 . AND. DIAPL.lt. 15) GO TO 20 6450 

IF(nPR0p2.NE.1) GO TO 25 6460 

IF(NCW.EO.I) YpRoP=TY(1.10) 6470 

IF(YY1 ( 1 ) .lE.YPROP.AND.YY(I.I) .GT.YPROP) L00P=L00P+1 6480 

IF (LUOP.lt. 2) GO TO 25 6490 

IF(NPRINT.EO.I.AND.NDIST.EQ.O) WRlTE(6.30) 6500 

30 FORMmKii*** TRAJECTORY TERMINATED. PROPELLER SLIPSTREAM "» 6510 

«, "ENTRAINMENT ***") 6520 

GO TO 20 6530 

25 CONfiNUF 6540 

IF(YY(1,3) .GT.YYl (3) )G0 TO 12 6550 

GO TO U 6560 

12 IF(ABS(yy(1.1 > ) .GT.1.EO)GO TO 13 6570 

GO TO l4 6580 

13 IF(YY(1 ,3) .6T.ZZ0.AND.ABS (YY(1, 1 ) ) .LT.ABS(YY1 (I) ) ) GO TO Id 4590 

14 CONTINUE 6600 
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6610 

WRIlfc. (6,15) 6620 

15 FORMAT(.IO"»>'»**« DO lOUP PARAMETER EXCEEDED IM SUB2D «***") 6630 

GO TO 2n 6640 

16 CONTINUE 6650 

iF(NUIST.EO.O.AMn.NPRlNT.eo.l) !^RITE (6.17) f.660 

17 FORMAT('inii,ii«**o trajectory COMPLETED. Z LESS THAN ZERO **»*••) 6670 

call YZFIN(Yy.YYl.NVAR.T.Tl.L) 6680 

GO TO 2o 66R0 

le COnTINUF 6700 

iF(NOIST.EO.O.ANn.NPRlNT.eQ.l) WOlTE (6,19) 6710 

19 FORMAT (mom, (,,»*, trajectory TERMINATED. VORTEX ENTRAINMENT *«*on) 6730 

20 CONTINUE 6730 

JJ=JJ+K 6740 

IRUN( IC0UNT)=K 6750 

Y=YY(1,i) 

V0=YY(1,5) 6770 

Z=YT (1.3) 6780 

wD=YY(1,4) 6790 

iF(NtVAp.EQ.l) OIA=YYd.NVAR)/l.E-6 6800 

iF(NtVAp.EO.O) DIA=0IA1/1.E-6 6810 

return 6820 

END 6830 

6840 

SUHROuTINF sUP3D; called 3Y the main PROGRAM. SU83D CONTROLS 6850 

The ToAJECTORY CALCULATION FOR ALL 3D CASES. GIVEN THE 6860 

initial VALU'^S and user selected OPTIONS. SUB3D CALCULATES 6870 

]HE TRAJECTORY AND REjuRNS THE FINAL VALUES OF THE VARIABLES 6880 

6890 

subroutine SU03D(NVAR.X.Y.Z.UO.VO.WD.T,OIA.YY,SAVE.CSAVE. 6900 

X YMAX.ERrOR.Oy.YYI.PW. IP) 6910 

dimension YY(13.NVAR) tS^VE (13.NVAR) .CSAVE (NVAR . 3 ) .YMAX(NVAR) 6920 

dimension error (NVAR) »DY (NVAR) ,yy1 (NVAR) 6930 

dimension PW(NVAR.NVAR) . ip (NVAR) , IRUN(50) 6940 

COMMON /aREAI/ A,CL.0»U,G*ZYO,7ZO.OA.VIS.VCW.D.DD 6950 

COMMON /AREA2/ N20 , N30 1 Ng y AP . FPS .NPR0P2. New . NTUN . ND I ST , NPR INT 6960 

COMMON /areas/ ZPROPa’P*^! .PK2.PK3.PK4.RP 6970 

COMMON /AREA6/ J J , I RUN . ICOUNT 6980 

COMMON /APEAP/ DIAl 6990 

7000 

INITtaLIZE the dependent variables and store in array 7010 

7Y(1,nVAP) 7020 

7030 

0IA=DIA*i,E-6 7040 

DIAI=0IA 7050 

C=3.14lsq2653589793 7060 

T=0.E0 7070 

yY(1,1)=v 7080 

YY(I.2)=vO 7090 

YY(1,3)=7 7100 

YY(I.4)=hd 7110 

YY(I,5)=x 7120 

YY(1.6)=ijD 7130 

IF(NEVAR.EQ.I) YY(1,NVAR)=DIA 7140 

IF (N30.F0.2) GO TO 4 7150 
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IF (NVAR.6E.12) GO TO 1 7160 

if(nvar.ge.io) go to a 7170 

lF(NHROp 2 .eo.n YY(1»7)=ZpR0P;> 7180 

GO ro 4 7190 

1 YY(i*12)3ZPROP2 7200 

YY(1»11)=0. 7210 

2 YY(lil0)aZZ0 7220 

YY(l*9)a-ZY0 7230 

3 YY(ltfl)=ZZ0 7240 

YY(1»7)=7Y0 7250 

4 CONTINUE 7260 

7270 

initialize values for DIFSuB 7280 

7290 

H=l.fc-4 7300 

MF=a 7310 

DO 5 I=1,NVAP 7320 

YMAX(I»=i.E0 7330 

5 CONTINUE 7340 

MAxDtR=6 7350 

HMIN=1.f-16 7360 

HMAX=l.E0 7370 

jST«HT=o 7380 

LOOP=0 7390 

YPfiOP=0. 7400 

7410 

perform the step integration by successive calls of DIFSUB 7420 

7430 

UlFStlR IS A FORTUOI SUBROUTINE TO SOLVE A SYSTEM OF N SIMULTANE 7440 

UIFFEPENTIAL EQUATIONS 7450 

7460 

DO 1** K=1,1000 7470 

DO 9 I=1,NVAP 7480 

YYl (n=YY(l,I) 7490 

9 CONtINUE 7500 

Tl=i 7510 

7520 

CALL DlFSBM(NVAR.TtYY*SAVE.CSAvEtH,HMIN.HMAX* EPS «MF.YMAX .ERROR. 7530 

X KFLAG.JSTART.mAXDeR.PW. IP) 7540 

7550 

monitor variables after each TIME-STEP TO CONTROL OUTPUT 7560 

AND STORAGE OF TRAUEdORY AND TERMINATE TRAJECTORY IF 7570 

RARTicLE becomes entrained in the vortex OR HITS THE GROUND 7580 

7590 

0IAPL=Yy(1.NVAR)/1.E-6 7600 

iF(NtVAo.EO.l) YY(1 .NVAR)=yY(1,NVAR)/1.E-6 7610 

IF(nPRInT.EO.I.AND.NDIST.eq.O) write (6.11) K.T. (YY(l.I) .I=1.NVAP) 7620 
11 F0RMAT(.< ••,I4.4X.14F9.5) 7630 

IF(NtVAp.EQ.l) YY(1 .NVAR)=yY(1,NVAR)*1.E-6 7640 

L=K*JJ 7650 

IF(YY(1,3).lt.O.EO)GO to 16 7460 

IF(NUISt.EO.I) go to 10 7670 

WRITE(8,22) T.YY(l.l) ♦TY(1,3) 7680 

22 FORMAT (3F13. 6) 7690 

10 continue 7700 
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IF(NEVAP.EQ.i. and. DIAPU.lt. 15) GO TO 20 7710 

IF(NPROpp.NF.I) go to 25 7730 

IFINCW.fq. 1) YPR0P=YY(ltl0) 773O 

II- <YY1 (1) .lE.YPROP.AND.Y,',(1,1) .GT.YPROP) L00P=L00P + 1 7740 

IF(L00P.lT.2) go to 25 7750 

iFtNPRlMT.eo.l.AND.NOlST.EQ.O) WRITE(6.30) 7760 

30 FORMATCn*** TRAJECTORY TERMINATED* PROPELLER SLIPSTREAM '•* 7770 

i "FNTPAINMENT »**“) 77S0 

GO TO 2o 7790 

25 continue 7800 

IF (YY (1 ,3) .gt.YYI (3) )(jO to 12 7810 

GO TO I4 7820 

12 lF(Abs(YY(l.l) ) .GT.1.E0 )Go TO 13 7830 

60 TO l4 7940 

13 IF(YY(1.3) ,GT.ZZO.AND.A0b(YY(l,l)) .LT.ABSCYYI (in > GO TO 18 7850 

14 CONjlNUp 7860 

K=K“1 7870 

WRIIE (6.15) 7380 

15 FORMAT(iiO"*"»«** do loop parameter exceeded in SUB30 ****") 7890 

GO TO 20 7900 

16 continue 7910 

IF(NUISt.EO.O.ANO.NPRINT.EQ.I) write (6.17) 7920 

17 format (mom. 'H**** trajectory COMPLETED* Z LESS THAN ZFRO ****•>) 7930 

call YZFIN(YY*YY1*NVAR*T*T1*L) 7940 

GO TO 20 7950 

18 continue 7960 

IF(NUISt.EO.O.ANO.NPRINT.EQ.I) write (6*19) 7970 

19 FORMAT ("O"*'***** TRAJECTORY TFRHINATED* VORTEX ENTRAINMENT ****») 7980 

20 CONllNUE 7990 

JJ=JJ+K 8000 

IRUN(IC0UNT)=K 8010 

Y=YY(l,i) 9020 

VO=YY(l,p) 3030 

Z=YT(1*3) 9040 

wO=YY(1.4) 8050 

X=YY(1,5) 9060 

UD=YY(1,6) 8070 

IF(NEVAd.EO. 1 ) DIA=YY(1*NVAR)/1.E-6 8080 

IF(NLVAP.EQ.O) 0IA=0IA1/1.E-6 8090 

RETURN 8100 

END 8110 

8120 

PLOErv is required by OIFSUB. but not needed for this method 8130 

8140 

SUBMOUTInE PEDERV(A.0*C*O) 8150 

return 8160 

end 8170 

8180 

DIFFun is a SUBROUTiNt CALLED BY OIFSUB TO EVALUATE THE hICjHT- 8190 

HAND fide OF THE SYSTEM OF DIFFERENTIAL EQUATIONS 3200 

RHS IS RETURNED TO DlFsua IN ARRAY DY(NVAR) 3210 

8220 

subroutine DIFFUN(NVAR*T*YY*DY) 8230 

DIMENSION YY(13.NVAR) *0,T (NVAR) 8240 

COMMON /AREAl/ A*CL*B*U*G*ZY0,ZZ0*DA*VIS*VCW*0*DD 8250 
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COMMON /«REA2/ N2D*N3D»NEVAP*EPS»NPR0P2»NCW,NTUN«NDIST*NPHINT 326(1 

C0M_MO(\j /flREA3/ DSW3*PA»pV,SCN,0V 8270 

common /areas/ ZRR0P2»PK1 fPK2,PK3»PK4,RP 8230 

COMMON /AREA3/ niAl 8290 

lE(NevAP.EQ.O) DIA=0IA1 8300 

Y=Yt(l,i) 8310 

V0=YY(1,2) 8320 

Z=YY(1,3) 8330 

w0=YY(1,4) 8340 

IF(N30.F0.1) X=YY(1»5) 8350 

IF (NjO.ICO. 1 ) uD=YY(1»6) 8360 

8370 

OHOSSWIND MODEL 8380 

3390 

ZCw=A 8400 

IF (ZOw.lT.O. ) ZCW=0. 8410 

ZCi*=4CW*B 8420 

IF (NCW.CO.O) VCW=0. 8430 

8440 

HhANcH to correct SECTION OF CODE DEPENDING ON OPTIONS 3450 

IjLLEcTED 8460 

8470 

IF(N20.F0.2.0R.N3D.EQ«2) GO TO 45 8480 

iFINjO.Pfl. 1) GO TO 15 8490 

KK=0 8500 

IF (NVAR.FO.ll.OR.NVAR.EQ.lO) GO TO 10 8510 

IF(NL«i.f(3,1) go to 20 8520 

IF(NTUN.EO.l) GO TO 30 8530 

IF(NPr0p?.E0.1) go to 35 8540 

IF(NVAR.EQ.5.0R.NVAR.Em.4) go to 40 8550 

WRiTt (6.91 8560 

9 FORMAT (•• M,n***** NVAR INCORRECT *****'•) 8570 

15 CONTINUE 8580 

KK=<: 8590 

IF(NVAR.GE. 12) GO TO 10 8600 

IF (NVAR. gE .10) GO TO 20 8610 

IF(NVAR.gE.7.AND.NPR0P2.EQ.1) GO TO 35 8620 

IF (NVAR. GE. 6) 60 TO 40 8630 

WRITE(6,9) 8640 

10 CONTINUE 8650 

8660 

ROTh VORTICES and PROP POSITIONS ARE 8670 

UtTEPMINED BY STEP INTEGRATION 3680 

HtRE NCW=1 AND NPROP2=l 8690 

8700 

CALL VELTV(YY(1.5+KK) »TY(1,6+KK) ,YY(lt7+KK) »YY(lf8+KK) »Y,ZfVlfWl) 8710 

CAUL PR0P20(Y,Z.YY(1»'^*KK) ,yY(1.10+KK) fV2,W2) 8720 

VA=V1*V2*(VCW*ZCW«*.251/U 8730 

wA=wl*W2 8740 

IF(N30.eO.O) go TO 17 8750 

CALL 60uNO(X.Y,Z*UA.w3) 8760 

call PR0PX(Y.ZiYYn*Il) »YY(1*12) *UAP) 8770 

GA=UA*UAP 8780 

WA=>«A4W3 «790 

17 CONTINUE 8800 


77 



non ooo oooo 


CALL VElTV(YY(1,5+KK) »YY(1,6+kk).YY(1,7*KK) .YY(1.B+KK) . 881(1 

^ YY ( 1 «q+KK) . YY a . 10+KK) «VP1.WP1 ) 8820 

OY(Y»<K)=vP1+(VCW*(YY(1.10*KK)*8)**.25)/U 8830 

DY(1U+KK)=WP1 8840 

IF (YY(1,io+KK) .LT.RP) OY(10+KK)=0. 8880 

call VElC (YY ( 1 .5*KK) .YY (1 ,bi-KK) .YY<1 *7+KK) .YY (1 «8*KK) ♦ VltWl) 8860 

0Y(=*'<K)=V1*(VCW*(YY< l*6*KK)«P)**.25) /U 8 87 0 

0Y(O+KK)=Wl p0gO 

call VElC( YY< 1,7*KK) . YY (1 ,8*KK) ,yY( 1»5+KK) . YY (1 *6+KK) « V1*w1) 8890 

OY ( ^♦«K)=-vi* (VCW* (YY <l fO+KK)»0) **.25) /U 39 0 0 

DY(b+KK)=-Wl 8910 

GO TO 5n 8920 

20 continue 0930 

8940 

NC«=i! VOPTEX POSITIONS DETERMINED BY STEP INTEGRATION* 8950 

NO P0OP 8960 

8970 

CALL VElTV(YY(1,5*KK) ♦YY(1,6+KK) *YY(1,7+KK) *YY(I*8+KK) *Y*Z*V1*»/A) 8980 

VA=Y1+ ( vCW*ZCW**.25) /U 8990 

IF(N3D.fq.O) go to 25 9000 

CALL BOijnO(X*Y,Z*UA*w1) 9010 

«A=HA + ii/i 9020 

25 CONTINUE 9030 

CALL VElC (YY ( 1 *5+KK) * YY ( 1,6+Kk) , YY (1 *7+KK) , YY ( 1 *8+KK) . VI * >^1 ) 90 40 

DY(J+KK)=V1+(VCW*(YY<1*6*KK)*8)**.25)/U 9050 

DY(b+KK)=wl 9060 

CALL VElC(YY(1,7+KK) *YY (1,8+KK) *YY(1*5+KK) ,YY(1*6+KK) *V1.w1) 9070 

DY ( (+KK)=-Vl+( VCW* (YY( 1*8+KK)«0) **.25) /U 9080 

0Y(B+KK)=-W1 9090 

GO TO 5o 9100 

30 continue 9110 

9120 

ntun=i: nprop20=o or 1 9130 

9140 

VP1=0. 9150 

WP1=U. 9160 

CALL TUNVEL(Y,Z.YY(1 ,5) *YY(1 ,6) ,Vl*Wl) 9170 

IF(NHROP?.E0.1) CALL PROP20(Y,Z.O.*YY(1*7) ,VPl*<lPl) 9180 

VA=Vl+VPl+(vCW*ZCW«*.2b)/u 9190 

WA=wl*WPl 9200 

CALL V0pTUN(YY(l*5) ,YY(1*6) *0 y(5) *0Y(6) ) 9210 

IF(NPRoP2.EQ.O) go to 50 9220 

call TUNVEL(0. *YY(1 *7) «YY (1*5) ,YY(1*6> «VDUM,DY (7) ) 9230 

GO TO 50 9240 

35 continue 9250 

9260 

NHROp2=1s without CROsSWINO or tunnel 9270 

9260 

call VEL?0(T*0,.YY(1.5+Kk) *V1,W1) 9290 

IF (YY ( 1 ,8+kk) .LE.RP) rtl=0. 9300 

OY(P*KK)=w1 9310 

CALL VEL20 (T.YtZfVA. wA) 9320 

call PRnP20(Y.Z.O. *YY(1*5*KK)*VPI*WP1) 9330 

VA=VA*VPl 9340 

WA=WA+<ipi 9350 
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IF (N3O.Ff5.0) GO TO 37 9360 

CALL G0UMn(XtY»7.UA*w2) 9370 

y CALL PROPX(Y*Z»0.*YY(1»7) ,UAP) 9380 

gAsUA+uap 9390 

V,a=«iA*w?. 9AOO 

27 CONTINUF 9410 

GO TO So 9420 

40 CONTINUF 9430 

9440 

N<;o=i OR n30=i: no Options fxcept nevap may equal l 9450 

9460 

CALL VEL20(TtY*Z.VA,wA) 9470 

IF(N30.fq.O) GO TO 50 9480 

CALL 50uND(X.YtZ,UA»wl ) 9490 

V,A=«A + wi 9500 

GO TO 5o 9510 

9520 

N,jO=? OR N3D=? 9530 

OtT VELOCITY FROM USER SUPPLIED SUBROUTINE USERV 9540 

9550 

45 continue 9560 

CALL USFRV(T»X.Y,ZtUA*VA,WA) 9570 

50 CONIINUf 9580 

9590 

lalculate rhs of droplet Dynamics equations and evaporation 96oo 

9610 

IF(NLVAP,EQ.1) ni A=YY < 1 *NVAR) 9620 

R=U* ' (0A*DIA)*S0RT( (VA-VD) **2+ (WA-*(0)**2) )/VIS 9630 

IF(N30.EQ.1) R=u*{ <0A*DIA)*SQRT ( (UA-UD)»*2*(VA-VD)*«2 9640 

h +(WA-wO)*“2) )/viS 9650 

CDR=1.0+0. 197*R»*0.63+0,26E-3*R»*1.38 9660 

S=lb,*(q«0A) /(DIA«0D) 9670 

RU=(DA*oiA*U)/VIS 9680 

IF(NEVAP.EQ.O) GO TO 60 9690 

0Y(NVAR)=-2.*(B/U)“( 16.016/28.97) *(DV/OIA)*(DA/DD)*( (PSW8-PV) 9700 

X /(PA-PV) )*(2.+.6*SCN**(1./3.>*R**.5) 9710 

60 continue 9720 

0Y(I)=Yy(1,2) 9730 

DY(2)=( {5*CnR)/RU)*<VA-VU) 9740 

DY(3)=YY(1.4) 9750 

. DY(A)=( (S«CDR)/RU)*(WA-wD)-(B«G)/U**2 9760 

iF(NtVAP.EQ.l) DY(2)=0Y(2)-<3./0IA)*VD*0Y(NVAR) 9770 

IF(NEVAp.EO.I) OY(4)=OY(4)-(3./oiA)*WD*OY(NVAH) 9780 

IF(N3D.EQ.0) go TO 70 9790 

DY(5)=Yy(1.6) 9800 

DY(0)=( (S*COR)/RU)*(UA-UO) 9810 

IF(NtVAp.EQ.l) nY(6)=DY(6)-(3./0IA)*U0*DY{NVAR) 9820 

70 RETURN 9330 

END 9840 

9850 

iURRouTiNE USERV USER SUPPLIED 9860 

9870 

WITH N?0=? OR N30=2 USERV IS CALLED BY DIFFUN FOR THE 9880 

VELOCITY IN THE WaKE. THIS ALLOWS THE USER TO REPLACE 9890 

THE CODE'S wake MODEL WITH ONE OF THEIR OWN. THE NEVAP 9900 
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option i<; still available as is NOIST, NPRINT, and NPLOT. qgio 

The POOPELLER AND TUNNel MODELS ARE AVAILABLE TO THE 9R?0 

USER BY INCLUDING TmE PROPER INPUTS AND SUBROUTINE CALLS. 9930 

9940 

PmSSFO TO USE°V IS the NONOIMNSIONAL TIME.Tt AND THE 9950 

NUNOimensIONAL position X«Y,Z. the user RETURNS THE 9960 

LORResPONOING NONOIMEnSIONal VELOCITIES UA*VA»wA. 9970 

9990 

SUSHOUTINE USERV(T.X,T*Z*UAfVA,WA) 9990 

RETUHN 10000 

END 10010 

10020 

10030 

!>UBRouTINE VELPD: input DIMENSIONLESS TIME AND POSI T ION ( Y i Z ) 10040 

VtL2D RETURNS 0 1 MENS I ONLESs INOUCtO VELOC I T lES ( U A , V A ) FROM 10050 

M VOrtEX PAIP unaffected By CROSSWINO ioobo 

10070 

SUBHOUTinE VFL2D (T.Y,Z*VAfWA) 10080 

COMMON /AREAI/ A,CL*B»U.U,ZY0.ZZ0»DA.VIS»VCW.DiD0 10090 

10100 

LALCuLATE the position OF the vortices 10110 

10120 

C=3* 1A1592654 10130 

C1=1/(Zy0**?)+1/(ZZ0*«2> 10140 

C2=-t (C1*ZY0**2-2.0)/SQHT(C1*ZY0**2-1.0J ) 10150 

C3= < (Ci«ZZO**2-2,0)/SQrT(C1*ZZO**2-1.0) ) 10160 

C4=C4*{(4.E0*C*A)/(CL*C1)) 10170 

C5=C3*( (4.E0*C*A)/(CL*C1) ) 10180 

IF(T.GT.C4) po TO 10 10190 

ZY=SURT( ( ( { ( (T<*CL*Cl)/t4.0*C»A) )»*2)*C1-(C2*T»CL*C1**2)/(2.0*C* 10200 

^ AJ+ci«C2*«2+4,0*C1)-Sqrt( ( ( ( (T*CL*C1)/ (4.0«C*A) )**2) *C1-(C2 10210 

H *I*Cl*C1**2) /(2.0*C*a)+C1«C2**2+A.o*C 1 )**2-A.0*Cl**2* ( ( (T*CL*C1 10220 

(, (A.o*C*A) )** 2-(C2*T«CL»Cl>/{2.0*C»A)+C2**2+4,0) ) )/(2.0 10230 

X *L1**2)) 10340 

GO TO 2n 10250 

10 zy=P0RT( ( ( ( ( (T*CL*C1)/(A.0*C«A) )**2>*C1-(C2*T*CL*C1**2) /(2.0*C* 10260 

X A) +Cl*C2««2+4.0*Cl ) *SURT{ ( t ( (T»CL*C1) /(4,0*C*A) )«*2) *C1-(C2 10270 

8, *T*Cl*C1**2) /(2.0*C»A)+C1*C2»*2+4.0*C1)**2-4.0»C1**2»( < (T*CL*C1 10280 

X l/(4.n*C*A) 1**2-(C2*T»CL<‘C1)/(2.0*C*A)+C2**2+4.0) ))/<2.0 10290 

X *L1**2)) 10300 

20 IF(T.LT.CS) GO TO 30 10310 

ZZ=SURT( ( ( ( ( (T*CL«C1)/<4.0*C»A) )**2)*C1-(C3*T*CL*C1**2)/(2.0*C* 10320 

X Aj +ci*C3*»2+4.0«C1»-SQRT( ( ( ( (T*CL*C 1 ) / ( 4 . 0*C*A ) )**2)*C1-(C3 10330 

X «T«Cl*CU*2) /(2.0»C*A) ♦C1*C3»*2+4.0*C1 )**2“4.0*C1**2*( ( (T*CL*C1 10340 

X J/ (4,fl*c*A) )»*2-(C3*T«CL*C1)/<2.0*C»A)+C3**2+4.0) ) )/ (2.0 10350 

X *L1*02)) 10360 

GO TO 40 10370 

30 ZZ=SQRT( ( ( ( ( (T«CL*C1)/(4.0*C*A) )**2)*C1-(C3*T*CL*C1**2)/(2.0*C* 10380 

X A) ♦Cl*C3**2+4.0»Cl )*Sqrt( ( ( ( (T»CL*C1)/(4.0*C*A) ) **2) *C1-(C3 10390 

X *f*CL»Cl*»2)/(2.0*C*A)*Cl*C3*»2-t'4.o*Cl)**2-4.0*Cl**2*(( (T*CL*C1 10400 

X >Y(4.o«C*A) )*»2-(C3»T<‘CL*C1)/(2.0*C*A)4C3**2+4,0) )>/(2.0 10410 

i *Cl**p)) 10420 

40 continue 10430 

10440 

LALCuLATE the induced velocities 10450 


80 



o o r> r> o n nonoo noon 


C 10460 

VA=(CL/(a.0*C*Al )*( uz*2)/l (Z7+Z)**2*(ZY-Y)**2)-(ZZ+Z)/( (ZZ+Z) 10470 

S. **2+(ZY*y)**3)-(ZZ-Z)/< (ZZ-Z)»»2+(ZY*Y)**2) +<ZZ-Z>/ ( (ZZ-Z)**2 10480 

& *(ZY-y)*«21) 10490 

WA=ICL/(2.0«C*A) )*( (ZY-Y)/( (Z7*Z) **2+ (ZY-Y) **2) + (ZY^Yl / ( tZZ*Z> 10500 

f, **2* (ZY*Y) **?)-( ZY-^Y) /< <ZZ-Z)**2*(ZY+Y)**2)-(ZY-Y)/( <ZZ-Z)»*2 10510 

i ♦<ZY-y)**2) ) 10520 

RETUHN 10530 

ENO 10540 

10550 

SUBROUTINE round: CALCULATES THE INDUCED VELOCITIES (UA.WA) 10560 

UUE TO THE BOUND VORTEx AT A POINT (XfY*Z> 10570 

10580 

subroutine ROUND (XtY.ZtUA.wA) 10590 

common /AREAl/ A.CL»0»U»GfZYO.7ZOfOA»VIS»VCWfO»OD 10600 

R5=SORT(x**a+(ZZO-Z)**2) 10610 

R6=5QRT(x**2+(ZZ0*Z)**2) 10620 

UA=l.-(CL/f4.*3,141592654*A) )*( ( (l.-Y)/SQRT(R5»*2-t-(l.-Y)**2)* 10630 

Si U.+Y)/S0RT(R5»*2+<1.+Y)**2) ) *( (ZZ0-Z)/R5**2>+( (l.-Y)/SQRT 10640 

t, «h6*«2+(l,-Y)<M>2) + (l. + Y)/SQRT(R6**2+a.'t-Y)**2) )*((ZZ0 + Z) 10650 

f, /H6*«2)) 10660 

WA=tLL/(4.*3.141592654*A))*(-( ( 1 ,-Y» /SORT t R5**2+ ( 1 ,-Y ) **2 ) ♦ 10670 

S ( 1 .■••Y) /SQPT (R5*«2* ( i .+Y)«*2) )» (X/R5»*2) + ( ( 1 .-Y) /SORT (R6**2 10 680 

S ♦« 1.-Y)**2)+(1.*Y)/SQRT (R6**2+(1.*Y)**2) )«(X/R6**2) ) 10690 

RETURN 10700 

ENO 10710 

10720 

SUBROUTINE PROPX: RETURNS THE X VELOCITY IN THE PROPELLER 10730 

SLIPSTREAM UAP DUE TO a PRQP CENTERED AT (YP*ZP) AT A 10740 

POINT (Y.Z) 10750 

10760 

subroutine PROPX (Y,Z»YP*ZP,UAP) 10770 

COMMON /AREAl/ A,CL»BiU»G.ZYO,ZZOtOA«VIS*VCW.O*DD 10780 

COMMON /areas/ ZPROP2*PKf *PK2,PK3*PK4,RP 10790 

Pl=buRT ( (Y-YP) **2+ (Z-ZP) «*2) 10800 

IE(RI.LT.RP) go to 10 10810 

UAP=0. 10820 

RETURN 10830 

10 CONTiNUr 10840 

R2=Rl/Rp 10850 

L'AP=(PKi*(1.-r2)+PK2*R2*5IN(R2*3. 141592654) )/U 10 860 

RETURN 10870 

END 10880 

10890 
10900 

subroutine PR0P20: returns OImENSIONLESS propeller INDUCED 10910 

VELOCITIES (V»W) AT POSITION (YtZ) FOR A PROP AT (YP«ZP) 10920 

(ALL POSITIONS DIMENSIONLESS) 10930 

10940 

subroutine PROP2n <Y»Z«TP.ZP»V, w) 10950 

COMMON /AREAl/ A,CL*8«U»6.ZYO,ZZO»OA.VIS»VCW»O.DD 10960 

COMMON /areas/ ZPR0P2*PK1»PK2,PK3*PK4,RP 10970 

Rl=b(JRT ( (Y-YD) **2+ (Z-ZP) «*2) 1098 0 

IF(RI.LT.RP) GO TO 10 10990 

V=0. 11000 
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w=o* linio 

RETUKN 11020 

10 CONTINUf 11030 

R2=Hi/Ro 11040 

SVl=(PKT»<i,o-P^)-t.PK4«R2*SIN(R2»3. 141592654) ) /U 11050 

V=( U-Zp)/R1)*SV1 11060 

W=( (YP-y) /PI ) *SV1 11070 

RETUHN 11080 

end 11090 

11100 

11110 

Subroutine veltv: returns inOuCEd velocities (v,w> at iii?o 

POSITION (Y,z) FOR VORTICES LOCATED AT (ZYR*ZZR) AND 11130 

}ZYL.ZZL) (ALL QOANTITIES DIMENSIONLESS) 11140 

11150 

SU8«0UTinE veltv (ZYR,ZZR fZYL»Z7L*Y.Z.V.W) 11160 

COMMON /4REA1/ A,CL«0»UtG.ZYO,ZZO*OA»VIS» VCWtO.DD 11170 

COMPLEX 0C,ZCR.ZCL.WZ»ZD 11180 

G AMMA=Cl/(A*p. *3. 14 15*^2654 ) 1 1 lOO 

GC=CMPLX (0. .gamma) 11200 

ZO=OMPLx (Y.Z) 11210 

ZCR=LMPlx(ZYR.ZZR) 11220 

ZCL=CmPlx (Zyl«ZZL) 11230 

WZ=-OC* (1 ./(ZD-ZCR)-I ./ (ZD-ZCL)*1. /(ZD-CON JG(ZCL) ) 11240 

*. -i./(7D-C0NJG(ZCR) ) ) 11250 

V=R£AL(W7) 11260 

W=-A1MAG(WZ) 11270 

retorn 11280 

END 11290 

11300 

11310 

subroutine VELC; returns induced VELOCITIES (V.W) ON A 11320 

VORTEX CORE AT (Y1»Z1) AND (Y2»Z2) (ALL QUANTITIES ARE 11330 

UIMEnrIONlESS) 11340 

11350 

subroutine VELC(y1*Z1»Y2.Z2.V.W) 11360 

COMMON /aREAI/ A,CL»B«U,GtZYO,ZZO*OA.VIS» VCW«D»DD 11370 

COMPLEX GC.ZCl *ZC2»WZ 11380 

GAMMA=Cl/(2.*3.141592654*A) 11390 

GC=CMPLx (0.. gamma) 11400 

ZC2=CMPlx (Y2.Z2) 11410 

ZC1=CMPlx(Y1,Z1) 11420 

WZ=-GC»(_1./(ZC1-ZC2)*1./(ZC1-C0NJ6(ZC2) ) -1 . / ( ZC 1-CONJG ( ZC 1 > ) ) 11430 

V=REAL(wz) 11440 

W=-A1MAg (WZ) 11450 

return 11460 

END 11470 

11480 

11490 

SUBROUTINE VORTUN: returns INDUCED VELOCITIES (VV.xv) 11500 

AT A point (ZYtZZ) FOR THE TUNNEL CASE WHERE VORTICES 11510 

ARE located AT (ZY»ZZ) AND (ZY»-ZZ) (ALL QUANTITIES ARE 11520 

DIMENSIONLESS) 11530 

1 1540 

subroutine VORTUN(ZY,ZZ. VV.WW) 11550 


82 



uuu uuooo uou 


COMMON /ARE41/ A,CL»0*U*G.ZYO,ZZOtD4.VIS*VCW»O*DD 11560 

complex Zl.w , const 11570 

C=3*i'^l5<52654 11580 

DTUN=D/R 11590 

Zl=CMPLx (ZY,7Z) 11600 

GAMMA=Cl/(2.*C*A) 11510 

E=0. 11620 

CONS>T=CmpLX (EtGAMMA) 11630 

W=(0..0.) 11640 

11650 

calculating the COMPLEX velocity at a vortex 11660 

11670 

00 10 1=1.21 11680 

N=I-ll 11690 

IF(N.EQ.O) GO TO 10 11700 

W=W*LONst*(-(1./(N«OTUN) )+(1./(2.*ZY+N*DTUN) )-(1./(2,*Z1 11710 

& *N*DTUN))+(U/(Z1-C0NJG(Z1)+N*0TUN))) 11720 

10 CONTINUr 11730 

W=W+CONst*( n./(2.*ZY) )-(l./(2.*Zl) )Ml./(Zl-C0NJ6(Zin ) ) 11740 

VV=MtAL(W) 11750 

HW=-AIMaG(W) 11760 

RETU«N 11770 

ENO 11780 

11790 

TUNVEl solves for the induced velocity (VA.KA) at a POINT IN TH 11800 

tUNNEL (YO.ZO) and a given VORTEX POSITION (ZY.+-ZZ). 11810 

‘All quantities are dimensionless) iiezo 

11830 

SUBHOUTimE TUNVEL (YO.ZD.Zy.ZZ.VA.WA) 11840 

COMMON /AREAI/ A.CL.8*U.G,ZYO,ZZO.DA.VIS»VCW.O»DO 11850 

COMPLEX C.Z.CONST.W 11860 

cc=3. 141592654 11870 

GAMMA=(cl*U*B)/A 11880 

C0=6AMMa/(2.«CC) 11890 

E=0. 11900 

CONST=CmplX (E.CO) 11910 

F=CC/0 11920 

OYO=YO*S 11930 

DZD=ZD*R 11940 

0ZY=ZY*8 11950 

OZZ=ZZ*R 11960 

C=CMPLX (OZY.OZZ) 11970 

Z=CMPLX (OYO.DZD) 11980 

W=-CONSt*( ( fCSIN(F*Z) )**2-(CSIN(F*C0NJG(C) ) )**2)/ ( (CSIN(F*Z) )**2 11990 

t. -<CSIn(F*C) )**2) )*( ‘ (2.*F*CSIN(F*Z)*CC0S(F*Zn*( (CSIN(F*Cn**2 120 00 

-(LSlNfF*CONJG<C) n**2) ) / ( (CSIN(F*Z) ) **2- (CSIN <F«C0NJG <C) ) )**2 12010 

i )«»2) 12020 

W=CONJG(W)/U 12030 

VA=KtAL(w) 12040 

VlAsAlMAGlw) 12050 

RETUHN 12060 

END 12070 

12080 

subroutine swath: calculates the distribution from EACH nozzle 12090 

AND SUMS THESE UP TO GENERATE A TOTAL DISTRIBUTION AND 12100 
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ti>TlM4TES THE drift. IF NPRINT=1 EACH NOZZLE'S DISTRIBUTION 12110 

AND nPIFT IS OUTPUT 12120 

12130 

SUBHOUTinF SWATH<NCOL*NROw,Q,0IAMN«STOEV*YNOZ*DWIDTH. 121A0 

^ NRUNtMET. width, OIST. DRIFT) 12150 

CIMtNSION TRAJ( 15,100) , YNOZ(lOO) , SIZE ( 1 5 , 1 0 0 ) ,NRUN( 100) ,OIST(l) 12160 

dimension ORIFTnOO),YI(l5),C{4,l5),SIZ«15),SlZl(15) 12170 

DIMENSION 0IAG( 15,100) ,S ( 1S,5) ,MET(100) ,XLIST (A) ,SLIST(A) 12150 

OlMtNSION 0IAMN( 100) ,STOEV(10O) ,0(100) ,C1 (4, 15) ,DLIST (A) 12190 

external CINT,8INT 12200 

COMMON /4REA2/ N2D,N3D,NEVAP,EPS,NPR0P2,NCW,NTUN,NDIST,NPRINT 12210 

COMMON /4REAQ/' SIZ,SIZ1,NR,C1,0IAMNI,STDEVI 12220 

COMMuN /4REA10/ TRA J , 0 I AG , Si ZE 12230 

DATA S/75*0./ " 12240 

KJ=0 12250 

STAHTs-wiDTH 12260 

NWlUfH=ln0»WIDTH+l-START*100 12270 

DO 10 I=i,nCOL 12280 

NR'=NRUN(I) 12290 

IF(NR.EO.O) 00 TO 200 12300 

IF(MtT(l) .EQ.I.OR.NR.LT.4) go TO 100 12310 

12320 

IE Met(I)=0 THE I»TH NOZZLE'S DISTRIBUTION IS CALCULATED 12330 

USINf! A CUBIC spline 12340 

12350 

IF(NPRInt.EO.O) go TO 46 12360 

wRIJc. (6.40) YNOZ(I) 12370 

40 FORMAT (••1".«yN0Z="«F5.3) 12380 

WRITE (6,45) 12390 

45 format ( HQ", 3 ( 7X, "YG'*,1 IX, "D I A",5X, "CONCENTRATION") ) 12400 

46 CONTINUE 12410 

DO EU J=i,nR 12420 

K=NK-J+i _ 12430 

SIZlU)=srZE(K,I) 12440 

SIZI (J)=nlAG(K,I) 12450 

Yl (U)=l ./(TP4J(K,I)-YN0Z (I) ) 12460 

20 continue 12470 

ISC='/ 12480 

AA=U. 12490 

BB=U. 12500 

IF(NEYAp.EO.I) call CSPL(SIZ,SIZ1,NR,C1,S,IBC,AA,BB) 12510 

CALL C5PL(YI,SI7,NR,C,S,IaC,AA,88) 12520 

IF(AbS(TRAJ(NR,I) ) .LT.OWIDTH) 60 TO 21 12530 

DELTAY=i./{nwIDTH-YNOZ(I) ) 12540 

CALL CSfiNO(vI,SIZ,NR,C,OELTAY,ORIFT(I) ,DUMI,DUHY) 12550 

21 continue 12560 

12570 

calculate the DISTHI'^UTION and add TD OIST array, 12580 

Output to 6 if nprint=i 12590 

12600 

NPTS=(Tp 4J(1 , I)-START)*l00+l 12610 

NENU=NWtr)TH 12620 

IF(TH4J(NR.I).LT.TRAJ(1,I) ) NeNO=1 12630 

NSTAHT=nPTS 12640 

NST0H=NEND 12650 
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IF{NE,N0.mE.1) go to 25 12660 

NSTART=meMD 12670 

NSTOH=Ndts 12680 

d5 CONTINUE 12690 

DIAULO=inoon. 12700 

DO JO KkksnSTART.NSTOP 12710 

KK=KKK 12720 

IF(NENO.EQ.l) KK=NST0P-(KKK-1) 12730 

Y6=bTAHT't-(KK-l)*.01 12740 

OeLTAY=l./(YG-YNOZ(I) ) 12750 

CAUL CSFJNO (YI tSlZfNRtCtOELTAY.DlAtOOYG»DUMY) 12760 

IF(L'IA.gE.OIAOLD) WRITE(6»27) I,0IA 12770 

27 FORM«T{iion,ti»«* CUBIC SPLINE PROBLEM NOZZLE NO.»*I3»" OlA="» 12780 

6 E13.6." ***") 12790 

OIAOLD=OIA 12800 

IF (UIA.lT.ORIFT(I) ) GO TO 35 12810 

IF(U1A.gT.SIZ(NR) ) GO TO 30 ' 12820 

0V0D=(EXP (-.5*( (DIA-0IAMN{ I) ) /STOEV ( I ) ) **2) ) 12830 

(. /(STPEV(I)*S0RT<2. »3. 141592654) ) 12840 

IF(NtVAO.EQ.l) CALL CSF INQ ( SIZ tSI Zl tNR ,Cl tOI A *DI A 1 .DUMY »D0MMY ) 12850 

IF (NtVAp.EO.l ) OVOO=OVUd*(OIAi*»3/OIA**3» 12860 

DY6=-1./((YG-YN0Z(I) )**<^) 12870 

DOYG=DDyG»OYG 12880 

OIST1=Ar<;(Q(I)«OvDD*DOYG) 12890 

0IST(KK)=0IST1+DIST(KK) 12900 

IF<NPHINT.EO.O) go to 30 12910 

IF(u1STi.LT.(1.E-4/Q(I)>) go TO 30 12920 

IF(KJ.E0,3) KJ=0 12930 

KJ=KJ+1 12940 

XLIST(KJ)=YG 12950 

0LIST(KJ)=0IST1 12960 

SLIST(Kj)=OIA 12970 

IF(KJ.E0.3) LRITE(6f80) < XL I ST (K 1 ) «SL I ST IK 1 ) . DL 1ST (K 1 ) »K 1 = 1 , 3 ) 12980 

80 FORMATI" ".3(E15.6t2El3.6) ) 12990 

30 CONTINUE 13000 

35 CONTINUE 13010 

lF(KJ.NE.3.AN0.NPRINT.tQ.l) WrITE< 6«80) { XL I ST (K1 ) t SL I ST ( K 1 ) , 13020 

S. DLIST(KI) «K1=1.KJ) 13030 

GO TO 200 13040 

100 continue 13050 

13060 

IF MeT(I)=1 the I'TH NOZZLE'S DISTRIBUTION IS CALCULATED 13070 

USING A linear fit 13080 

13090 

KJ=0 13100 

IFINR.LT.2) GO TO 200 13110 

IFINPRInt.EO.I) WRITE(6«10S) YNOZII) 13120 

105 FORM AT (•' 1 ", II YNOZ=" » F5 • 3/" 0",A(8X*'*YG*'»7X* "CONCENTRATION") /) 13130 

00 110 Jsl.NR 13140 

K=nR-J*i 13150 

SlZI'J)=SIZE(K,I) 13160 

ctZl (J)sOIAr,(K,n . 13170 

Yl I'J)=TRAJ IK, I ) 13180 

110 CONTiNUc 13190 

IF (NtVAP.EQ.O) GO TO 115 13200 


85 



noon 


lF(NR.Gfr.4) GO TO 112 13210 

WHlTt(6,ll3) I 13220 

113 FORMAT(ho"«i<*»* insufficient data to calculate distribution ***«/ 13230 

R, II *«* nozzle nO* "*I2.h distribution set = TO ZERO ***") 13240 

GO p 3oo 13250 

112 continue 13260 

I0C=0, 13270 

AA=U. 13280 

PBaO. 13290 

CALC CSpl (SIZ.SIZI «NR fCl ,StIBC.AAt88) 13300 

115 continue 13310 

13320 

calculate distribution and add to OIST array. 13330 

IT NPrINT=1 output to 6 13340 

13350 

YMIN=100 13360 

YMAX=-Inn 13370 

DO 1«:0 K=1.NP 13380 

IF ( Y1 (K) .LT. YHIN) YMIN=YI(K) 13390 

IF(Y1 (K) .GT.YMAX) YMAXsYkK) 13400 

120 CONTINUE 13410 

ISTART=(YMIN-START)*100+1 13420 

DO iJO K=ISTART.NWIDTH 13430 

YG=(K-1)«.01*START 13440 

IF(ABS(YG) .GT.DWIOTH) 60 TO 130 13450 

0ISTI=0. 13460 

IF ( YG.GT.YMAX) go to 199 13470 

J1=NK-1 13480 

DO 140 J=1,J1 13490 

IFIYG.LT.YI (J) .AND. YG.lt. YI (J+1) ) GO TO 140 13500 

IFIYO.Gt.YI (J) .AND.YG.GT.yI tJ+1) ) go to 140 13510 

SLOHE=(siZtj4l)-SIZ(jn/(YI (J+D-YI (J) ) 13520 

OIA=SIZ(J)+SLOPE*<YG-YI (J) » 13530 

DV0U=(EXP(-.5*( (DIA-DIAMN(I) l/STDEVd) )**2) ) 13540 

X /(ST0EV(I)*SQRT(2. *3. 141592654)) 13550 

IFINLVAP.EQ.O) GO TO 135 13560 

CALL CSfiND(SIZ.SIZ1.NH.C1,0IA,DIA1»0MY.0MMY> 13570 

DV0U=0VDD*(DIA1»*3/DIA**3) 13580 

135 CONTINUE 13590 

DISri=OlSTl+Q(I)*A8S(SLOPE*OVOO) 13600 

140 continue 13610 

DISTtK)=niSTl+DIST(K) 13620 

IF (NPRInT.EO.O) go to 130 13630 

IF(K0.E0.4) KJ=0 13640 

KJ=Kj+1 13650 

XLISUKJ)=YG 13660 

OLlSf (KJ)=DIST1 13670 

IF(KJ.E0.4) wrITE(6.141) (XLIST(Kl) .OLIST(Kl) .K1=1.4) 13680 

141 FORNATdi ".4 (E16.6.E14.6) ) 13690 

130 continue 13700 

199 continue 13710 

IE(KU.Ne.4.and.NPRINT.E0.U wRITE( 6.141) (XLIST(Kl) .DLIST(KI) 13720 

4 .K1=1.KJ) 13730 

200 continue 13740 

C ■ 13750 
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U«IFT ESTIMATIOM 13760 

13770 

ORFT=DRirT(I) 13780 

IF(0hifT(I),FQ.0) go to 300 13790 

IF(0RIFt(I).ME.SIZE(1’I> 1 GO TO 205 13800 

ORI>-T(I)-100. 13810 

GO TO 3so 13820 

205 continue 13830 

AA=sU. 13840 

0a=(DHlFT( D-OIamN (I ) ) /STOEV (I ) 13850 

ISEG=A8s (08*20) /2 13860 

lSEU=ISFr,*2 13870 

CALL SIMP ( AA.00. ISEG»CINT,ANS) 13880 

ANS=ANS*5QRT(1./(2.»3«1‘*1592654) ) 13890 

drift (I)=(.50+ANS)»100. 13900 

300 continue 13910 

IF(NtvAp.EQ.O) GO to 350 13920 

13930 

include the effect of evaporation in the drift 13940 

13950 

IF(NH.GE.4) go to 320 13960 

WRIIt(6,330) I 13970 

330 FORMAT («n"*'****INSUFFlCIENT DATA TO ESTIMATE DRIFT ***■■/ 13980 

(, M **» nozzle no. »»ia,H drift set equal to zero »**") 13990 

DRITTdjsO. 14000 

GO TO 10 14010 

320 CONTINUp 14020 

AA=URFT 14030 

P0=i>IZ(NR) 14040 

ISEOsIFix ( (SIZ (MR)-SIZ< 1) ) /20)»20 14050 

lFdSEG.LT. 2) ISEG=2 14060 

DIAMNI=niAMN (I ) 14070 

STOtVI=STOEV (I ) 14080 

CALL SImp(AA,0R,ISEG,OINT,ANS) 14090 

EVAP=1._aNS 14100 

drift d)=EVAP*iOO.+DRlFT(I» 14110 

350 CONTInUF 14120 

IF(NPRINT.EQ.I) WRITE <6t60) DRIFT(I) 14130 

60 format (hO"»"DRIFT="iE1A. 6) 14140 

10 continue 14150 

return 14160 

END 14170 

14180 

51HP IS A SIMPSON'S RULE INTEGRATION ROUTINE 14190 

14200 

AA, 0 p are the limits of Integration 14210 

ISEG IS THE EVEN NUMBER OF SEGMENTS AA TO BB IS DIVIDED INTO 14220 
UlNT IS THE subroutine WHICH EVALUATES THE INTEGRAND 14230 

ANS RETURNS the ANSWER 14240 

14250 

subroutine SIMP(AA.0B'ISEGtCINT«ANS) 14260 

DELTAs(rp-aa)/ISEG 14270 

IS£G1=ISEG-1 14280 

SUM=U. 14290 

DO IV I=1.ISFG1 14300 
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P®c='*. U 3 in 

K=I/2 U320 

IF( (^*2) ,EQ.n FAC=2. 14330 

X=AA'*I<*nFLTA 14340 

CALC. CInT(X,Y) 143<;0 

SUM=SU‘J*v*FAC 14360 

10 continue 14370 

CALL CInt(Aa,y1) 14330 

CALL CInT(‘Tf.y2) 14390 

SUM=SUM*y1+Y? 14400 

ANS=10ElTA/3. ) *<5UM 14410 

return 14420 

14430 

14440 

lint evaluates integrand when NEVAP=0 144S0 

14460 

subroutine CINTCXfY) 14470 

Y=EXP (-,i,0*X<n*2) 14480 

RETURN 14490 

ENO 14500 

- 14510 

BINT evaluates integrand RHEN NEVAP=1 " 14520 

14530 

subroutine RINT(X.Y) 14540 

DIMENSION SIZ ( 15) tSIZl (15) ,C1 (4.15) 14550 

common /areas/ SIZ.SIZ1*NR,C1,DIAMNI.STDEVI 14560 

CALL CSFIND(SIZ.SIZ1»NR, Cl, X.DIAI.OUMY. DUMMY) 14570 

DIA=X 14580 

DV0U=(ExP(-.5*{ (0IA-0IAMNI)/STDEVI)**2) ) 14590 

S. / (STDFVIOSQRT (2.03.141592654)) 14600 

Y=(DIAl/niA)*«3oOVOD 14610 

RETURN 14620 

END 14630 

14640 

INCOn evaluates the initial velocity of the droplet 14650 

14660 

Tu Correct for the effect of the bound vortex the droplet is 1467o 

given an initial velocity equal to its terminal velocity at 14680 

that point in the FLUR field 14690 

14700 

subroutine INC0N(YY.0IA,NVAR) 14710 

dimension YY(13,NVAR) 14720 

COMMON /aREaI/ A,CL.a«U.G.ZY0,zZ0.DA.VIS. VCW.D.DD 14730 

common /AREA2/ N2D.N3D,NEVAP,EPS.NPR0P2.NCW.NTUN»NDIST.NPRINT 14740 

C=3. 141802654 14750 

ZY=ZY0 14760 

ZZ=ZZ0 14770 

Y=YY(l,l) 14780 

Z=YY(1.3) 14790 

IF(NTUN.EO.I) go TO 5 14800 

VA=(Cl/(2.0»C«A) )*( (Z2*Z)/( (ZZ*Z)**2*(ZY-Y)**2)-(ZZ*Z)/( (ZZ+Z) 14810 

«, *"24(7Y*Y)**2)-<ZZ“Z)/< (ZZ-Z)**2*(ZY4Y>**2)+(ZZ-Z)/( (ZZ-Z)**2 14820 

X +(ZY-y)**2)) 14830 

WA=(Cl/(2.0*C*A) )* ( (ZY-Y)/( (ZZ+Z)**2*<ZY-Y)**2)* (ZY4Y)/( UZ*Z) 14840 

^ **2*(7Y*Y)«*2)-(ZY+Y)/( (ZZ-Z)*»2+(ZY*Y)**2)-(ZY-Y)/( (ZZ-Z)**2 14850 
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h +<ZY-y)*»2)) 14860 

5 CONTINUF 14070 

IF(NTUN.PQ.I) CALL TUNVEL ( Y t Z , Z Y. ZZ . V A . WA ) 14000 

S=10.* <B*0A) /(0IA*00) 14090 

RU= (DA*ni A«U) /VT'5 14900 

V0=VA 14910 

YYU,2)=vA 14920 

WDL=WA 14930 

WDR=WA 1<»940 

no 10 J=l,100 14950 

WOR=*iOR + WA 14960 

q=U«A0S( ( (DA*0 IA)*(WA-w0H) )/VIS) 14970 

C0R=i .0*0. l97*R«i*O.6340.26E-O3«R*«1.38 14900 

fr=( (S*cnR)/RU)»(WA-w‘^R)-(a*G)/u *»2 14990 

IF (I-R.Gt. O.FO) Go TO 4 O 15000 

10 COnTINOF 15010 

WRiTt (6.20) 15020 

20 format (' in"."**»* WOR exceeded DO LOOP IN INCON 15030 

xRITt <(^.30) 15040 

30 formatoiv'*""**** right Endpoint in incon incorrect **«*••) 15050 

YY(1»4)=o.EO 15060 

go to 90 15070 

40 CONTINUF 15080 

15090 

SOLVE FOR THE VELOCITY 0Y THE METHOD OF BISECTION 15100 

15110 

DO SO 1=1.500 15120 

WOM=(WOl+WDR) /2.E0 15130 

RsABS ( ( (OA*DIA) * (WA- wOM) ) /VIS) 15140 

R=R*U 15150 

COR=i.Q*0.197*R*»0.63*0.26E-03*R**1.38 15160 

FH=( (S*cnR)/PU)«(WA-W0M)-(B*G)/U**2 15170 

FF=ABS(FM) 15180 

IF(FF.LT.EPS)G0 TO 70 15190 

IF(FM.Gt.O.EO) wDR=wDM 15200 

IF(TM.Lt.O.EO) W0L=w0M 15210 

50 CONTINUF 15220 

wRiTt (6,60) 15230 

60 F0RMAT("1","**«* do loop parameter in incon exceeded ****") 15240 

YY(1,4)=o.EO 15250 

GO TO 9n 15260 

70 YY(i,4)=wOM 15270 

IF(NOISt.EO.O) write (6.80) VO.WOM 15280 

00 F0RMAT('in". "INCON V ALOES" , T 25 , "VO=" ,E 13 . 6. T45 » "WD=" , E 1 3 . 6// ) 15290 

90 RETURN 15300 

END 15310 

15320 

yzfin evaluates the droplet location and velocity 15330 

WHEN IT intersects THE GROUND PLANE 15340 

VALUFS FOR THE NEXT TO LAST TIME-STEP ARE IN THE ARRAY 15350 

YYl and T1. THE VALUES FOR THE LAST STEP YY AND T ARE 15360 

RtPLACEO RY THE linearly INTERPOLATED FINAL VALUES 15370 

TO 5f returned 15380 

15390 

subroutine YZFIN (YY.YTI .NVAR.t.TI »LI 15400 
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OIMtNSIOM YY (13.NVAR) « YYl (NV4R) , IRUN (50) 15*10 

common /AREAI/ A.CL»atU»G,ZYO,ZZO«OA«VISf VCW*0«D0 15*20 

COMMON /AREA2/ N20 . N30 1 NE V AP , EPS t NP90P2 1 NCW t NTUN » NO I ST . NP« INT 15*30 

common /AREAfi/ JJt IRUN» ICOUNT 15**0 

CONV=(O.EO-YY1 (3) ) / ( YT ( 4»3)-Yy1 (3) ) 15*50 

00 10 I=liNVAR 15*60 

YY ( 1 * I) = (YY (1 1 n-YYl (I) )*C0NV*YY1 (I) 15*70 

10 CONTINUE 15*80 

YY(lt3)=0. 15*90 

T=( > i-Tl ) *C0NV)+T1 15500 

IF(nOISt.EO.I) go to 20 15510 

WRlTt(8,30) T.YY(1»1) ‘YYllta) 15520 

30 F0RMAt(3F13.6) 15530 

30 RETUmn 155*0 

END 15550 

C 15560 

C subroutine CSPL 15570 

C CSPL FITS A CUBIC SPLINE THROUGH THE M INPUT POINTS (X*Y). 155B0 

C THE equation for THE CUBIC WHEN X IS GREATER THAN X(I) ANO 15590 

C LESii than X(I*1) is: 15600 

C Y(X)=C(l,I)»(X{I+n-X)**3*C(2,I)*(X-X(I))**3*C(3»I)* 15610 

C (X(I*l)-X)+C(*»I)*lX-X(ln 15620 

c 15630 

c pARAMFTfpS: 15640 

c 15650 

c X AN ARRAY INPUT AND DIMENSIONED XtM) CONTAINING THE 15660 

C INDEPENDENT VARI«B|.E IN INCREASING ORDER 15670 

C Y AN ARRAY INPUT AND DIMENSIONED Y{M) CONTAINING THE 15680 

C OF.PENDENT variable 15690 

C M number OF (X.Y) INPUT 15700 

C C an AROAY output WITH DIMENSION C(4*M-1> CONTAINING 15710 

C THE DESIRED COEFFICIENTS 15720 

C S A WORK ARRAY DIMENSIONED S(M*5> 15730 

C IBC =0 SECOND DERIVATIVE OF Y WRT X SET EJUAL TO A AT 157*0 

C X(l) AND 8 at X(M) 15750 

C =1 FIRST DERIVATIVE OF Y WRT X SET EQUAL TO A AT 15760 

C X(l) ANO 0 at X(M) 15770 

c A.B INPUT value of First or second derivative depending istbo 

C ON IRC 15790 

c 15B00 

C programed By M. RRAGG AARL/OSU MARCH*19B0 15810 

c 15820 

subroutine CSPL(X.Y,M»C.S.I0C,A.B) 15830 

DIMENSION X(l) ,Y(1) «C<^»1) *S(M,5) 15840 

DO R K2=lt5 15850 

DO 1 Kl=l,M 15860 

1 S(K1«K2)=0. 15870 

2 continue 15880 

IF(IOC.NE.O) GO TO 10 15890 

S(l»F)=i. 15900 

15910 

S(1»A)=a 15920 

S(M**)=B 15930 

GO TO 2n 159*0 

10 continue 15950 
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S(1»2)=3. 15960 

s«l*3)=l. 15970 

S(1*'*)=6.*((Y(2)-Y(1))/((X(2)-x(1) »**2) )-(6.*A/(X (2)-X(1) ) ) 15980 

S(M»1)=1. 15990 

S(Mi2)=?. 16000 

; (M*4)=6.* ( (Y<M-1 )-Y(M) ) / ( (X (M)-X<M-1) )**2) )+(6.*B/(X(M)“X (M-1 ) ) ) 16010 

20 CONjlNUc 16020 

Ml=M-l 16030 

DO 90 t=2»Ml 16040 

DX=AlI*i)_x(i) 16050 

DXl = X(I)_)((i_i) 16060 

s(i*i)=nxi/nx 16070 

Sa*<i) = (2.*DX*2.»0Xl)/0X 16080 

Stl»9)=l. 16090 

S(I*‘*1=6.*(({Y(I*1)-Y(I)) /DX**2)-< (Y (I)-Y (1-1 ) ) / (DX*DXl ) ) ) 16100 

30 CONTINUF 16110 

16120 

USE THOMAS ALGORITHUM TO SOLVE TRI-OIAGONAL SYSTEM 16130 


16140 

16150 

16160 

16170 

16180 

16190 

16200 

16210 

16220 

16230 

16240 

16250 

16260 

16270 

16280 

16290 

16300 

16310 

16320 

16330 

16340 

16350 


SCFINO USES THE X.Y.M.C FROM A CALL TO CSPL AND RETURNS FOR 16360 

any INOEPENOENT VARIABLE XXI THE CUBIC SPLINE VALUE OF THE 16370 

dependent variable YY* T«E first DERIVATIVE YP* AND THE 16380 

SECOND derivative Y2P* 16390 

16400 

SUBKOUTTnE CSFIND (X*Y«M»C»XX*YY»YP*Y2P) 16410 

DIMENSION X(l) ,Y(1) *C(A«1) 16420 

Ml=M-l 16430 

IF (X (1 ) .OT.X (M) ) GO TO 20 16440 

00 10 I=l«Ml 16450 

lF(Xm ,lE.XX.AN0.X(I*1) .GT.XX) GO TO 40 16A60 

10 continue 16470 

IF(XX.LT,X(n ) 1=1 16480 

IF (XX.Oe.X (M) ) I=M1 16490 

RO TO 40 16500 


DO 40 I=?,H 

S(I*‘i)=S(I»2)-S(I.l)*S<I-1.3)/S(I-l.2) 

S(I*A)=«;(I,4)-S(Itl)«S(I-lf4)/s(i-l,2) 
40 CONTlNUr 

S(M.b)=S(M,4)/S(M,2) 

DO bU I=?.M 
J=M-I+1 

S(J»P)=(S<J.4)-S(Jf3)*S(J+l,5))/S(J*2) 
50 continue 

DO 90 I=l»Mi 
DX=X U + 1 )-X ( I) 

C(1»I)=S(I*5)/(6.*DX) 

C(2»1)=S(I+1.5)/(6.*DX) 

C(3*I)=(Y(I)/DX)-( (S(I»5)*DX) /6.) 
C(4*I) = (y(I + 1)/DX)-( (S<I-t-1.5)*0X)/6.) 
60 CONIINUE 
betOhn 
END 

SUBROUTINE CSFIND 
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20 CONTINUF 16510 

DO 30 I=liMi 16520 

IF (* ( I ) .(iE.XX .ANO.X ( I + l ) .UT.XX) GO TO 40 16530 

30 CONTINUF 16540 

IF(XX.Gt.X(1) ) 1=1 16550 

IP(XX.UT,X (M) ) i=Ml 16560 

40 CONTINUF 16570 

YY=V <1 *I )« (X ( I + l )-XX) **3+C (2«I) * (XX-X( I) ) «*3*C (3 * I ) * ( X ( I + 1 ) -XX ) 16550 

& *U (4, I ) * ( xX-X ( I ) ) 16590 

YP=-3.*C(l«I)*(X(I+l)-XX)**a*3,*c(2*I)*(XX-X(I) )*®2-C(3*I) 16600 

X ■*■^-(4,1) 16610 

Y2P=o.*f;(i,i)*(x(I + l)-XX)+6.*C(2^n*(XX-X(I) ) 166 20 

return 16630 

END 16640 

16650 

bURPOUTINF DLOTl 16660 

16670 

►"LOTi IS rALLED IF NP(_OT=l AND PLOTS PARTICLE THAJECTOHIES 16660 

If NniST=0 AND THE FINAL DISTRIBUTION IF NOIST=X 16690 

16700 

SUBHUUTtnF PLOTl (TITLE) 16710 

OIMC.NSION IRDN(50) .YPLOT(lOOl) ,ZPLOT(1001) »TITLE(13) 16720 

COMMON /AREA2/ N20»N3D*NEVAP»EPS*NPR0P2»NCW»NTUN*N0ISTtNPRlNT 16730 

COMMON /AREA6/ JJ , I RuN t I COUNT 16740 

NT0T=JJ 16750 

rewind R 16760 

CALL PLOTF(120.4) 16770 

IF(NUISt.E0.1 ) GO TO 200 16780 

DO 10 IsI.nTOT 16790 

REAU(8,16) T.YPLOT(I) *ZPL0T(I) 16800 

16 F0RMAT(3F13.6) 16810 

10 CONTINUE 16820 

ZMAX=o. 16830 

DO 40 1=1 ,NT0T 16840 

IF(4PL0T(I) .GT.ZMAX) 2MAX=ZPLOT ( I ) 16850 

20 CONTINUF 16860 

YPLOT(NTOT*1)=0. 16870 

NT0T1=NT0T*1 16380 

CALL SCAlE(VPL0T,8.«NT0T1,FY.0Y) 16890 

FZ=U. 16900 

DZ=U7 16910 

IZ=ZMAX/0Y*1. 16920 

IF(IZ.LE.6) go to 30 16930 

CALL scale (ZPL0T.6. *NTOT»FZ»OZ) 16940 

DY=UZ 16950 

30 continue 16960 

AIZ=IZ 16970 

call PLOT(0.5tl.5«-3) 16980 

CALL AXIS (0. ,0. ♦ <2Y/B» *-4,8., 0. .FY»DY) 16990 

CALL AXIS (0. ,0., »2Z/B* »4 *AIZ»90..FZ»OZ) 17000 

CALL SYMROL(1.0«6.5f .1A«TITLE,0,«36) 17010 

K=1 17020 

00 ‘fU I=l.ICnuNT 17030 

II=IRUN( I ) +K-1 17040 

DO PO J=K«II 17050 
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JJ = J-K + 1 17(1^0 

YPLUT (JJ)=YPL0T (J) 17070 

ZPLOT(JJ)=ZPLOT(J) 170«0 

50 CONTIN'jr 17090 

NLAbT = lDijM<T) 17100 

call LlNP(YPLOTtFYfOY»ZPLOT.FZ.OZ*NLAST*0.1) 17110 

K=K*1PUn(I) 17120 

40 CONT l.NUF ' 17120 

CALL PLCTdl., -1,5*999) 17140 

RETUhn 17150 

200 CONTIMUF 17150 

YMAX=-10o. 17170 

YMIN=10n. 171P0 

DO ^;<;0 1 = 1, JJ 17190 

REAU18*?10) YPLOT(I) *ZPL0T(I) 17200 

2J0 FORMAT (2F13. 5) 17210 

IF (ZMLOt { I) .GT. l ,E-2 .ANO*YPLOt (I) *LT.Y“IN) YMIN=YpLOT ( I ) 17220 

IF (ZPLOT ( I ) ,GT. 1 ,E-2. AnO. YPLOT ( I) .GT.YMAX) YMAX=YPL0T ( I ) 17230 

220 continue 17240 

call PLnT(0.,1.5,-3) 17250 

K=0 17260 

DO 410 1 = 1 *j,j 17270 

IF(YPLOt{I) .LT.YMIN) go to 210 172«0 

IF ( YPLOt (I ) .GT.YMAX) OO TO 2l0 17290 

K=K+i 17300 

YPLOT(K)=YPLOT(I) 17310 

70LUT(K)=ZPLnT(I) 17320 

210 continue 17330 

call scale (YPLOT. a.*K»FX,OXJ 17340 

call SCalE(ZPL0T.6.,K»FQ,DD) 17350 

CALL AXl<; (0.,0. . 129/8' *-4,8., 0.*FX, OX) 17360 

call 4XTS(0.,0.,*CONCtNTRATION*,13,6.,90.,FO,OD) 17370 

call SYmPOL(1.0,6.5, ,U,tITLE,0.,36) 17380 

CALL LiNP (YPLOT,FX,DX*ZPLaT,Fn,DD,K,0,l) 17390 

CALL PL0T(11 .,-1.5,999) 17‘*00 

RETURN 17410 

END 17420 
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C-» This SUbtJOuTINe INTSGPATeS A SET OF N ORDINARY DIFFERENTIAL FIRST * 
C* ORDER equations OVER ONE STEP OF LENGTH H AT EACH CALL. H CAN SE * 

c* specified ry thf user for Each step, but it may be increased or * 
c* decreased py oifsur within The range HMIN to hmax in order to * 
c* achieve as large a step as possible while not committing a single * 

C« STEP ERROR WHICH IS LARGER IHAN EPS IN THE L-2 NORM. WHERE EACH * 
C» COMPONENT OF THF error IS DIVIDED BY THE COMPONENTS OF YMAX. * 
C» * 
C* THE PROGPAm REOUIRES FOUR SUBROUTINES N«MED * 
C* DlFhUN(N.T.Y.DY) --USEH-SuPpl lEO — 

C* OECUHP(N,M.pw,lP) — renamed to NOIOlZ IN THIS SOURCE — * 
C* SOLVE(n.m.PW,CSAvE(1.1> »IP) — RENAMED TO NOI02Z — » 
C* PEOEHV (T. Y.PW.M) --user-supplied— * 
C* THE FIRST. OIFFIJN. EVALUATES THE DERIVATIVES OF THE DEPENDENT * 
C* variables STORPD IN Y(l.I) FOR I = I TO N. AND STORES THE * 
C» OERIVAIIVES IN THE ARRAY DY. tHF NEXT TWO ARE CALLED ONLY IF THE * 
C* METHOD flag mf IS SET TO I OR 2 FOR STIFF METHODS. DECOMP IS A * 
C* STANOAhO LiJ decomposition with PIVOTING THAT DECOMPOSES THE MATRIX « 
C« PW, leaving thf PIVOTS IN THE INTEGER ARRAY IP. M IS THE DECLARED * 
C» SIZE Of Pw. IP(N) IS SET TO 0 IF PW IS SINGULAR. SOLVE PERFORMS * 
C* BACK substitution ON THE CONTENTS OF CSAVE(I.l). LEAVING THE * 
C* RESULTS THFdE. 

C* PEDERV IS USED ONLY IF mF Is 1. AND COMPUTES THE PARTIAL * 
C» DERiVAlIVEs OF THE DIFFERENTIAL EQUATIONS AS DESCRIBED UNDER THE « 

C« mF parawetfr. « 

C» * 
C» THF temporary storage SPACE IS PROVIDED BY THE CALLER IN THE * 
C» INTEGER array IP. THE ARRAY PW, AND THE * 
C* ARRAYS SAVE AND CSAVE. THE ARRAY PW IS USED ONLY TO HOLO » 
C* THE MATRIX OF THE SAME NAME. AND SAVE IS USED TO SAVE THE VALUES * 
C» OF Y IN CASE A STEP HAS TO BE REPEATED, BUT CSAVE IS USED TO HOLD * 
C* SEVERAL arrays. THE REGIONS USED ARE « 
C* CSAVt(I.i) IS USED mainly TO HOLO THE CORRECTION TERMS IN THE » 
C* corrector LOOP, and HOLDS THE DERIVATIVES DURING « 

c« JACOBIAN Evaluations. « 

C* CSAVt(I.p) IS USED TO SAVE THE VALUES OF THE SUMS oF ALL OF THE* 
C* CORRECTION TERMS IN THE PREVIOUS STEP AFTER THEY * 
C* have BEEN ACCUMULATED IN THE ARRAY ERROR IN THE * 
C* CURRENT STEP. THIS ENABLES THE BACKWARDS DIFFERENCE* 
C* OF error to be formed. IT IS USED TO ESTIMATE THE * 
C* STEP SIZE for one OROER HIGHER THAN CURRENT. « 
C* CSAVt(I,3) IS USED TO STORE THE DERIVATIVES WHEN THEY ARE » 
C» COMPUTED by OIFFUN. * 
C* * 
C* THE PARmmEtpRS to THE SUBROUTINE 0IFS8M HAVE * 
C* THE FOLLOWI^JG M^a^JI^4Gs.. * 
C* » 
C* N THE NUMBER OF FIRST OROER DIFFERENTIAL EQUATIONS. N * 
c* may he decreased on later Calls if the number of * 
C* ACTIVE equations REDUCES. BUT IT MUST NOT BE * 
C* INCREASED WITHOUT CALLING WlTH JSTART =0, * 
C* T THE independent VARIABLE. * 
C* Y A 13 RY N array CONTAINING THE DEPENDENT VARIABLES AND * 
C* their SCALED DERIVATIVES. Y(J+1,I) CONTAINS * 
c* the j-th derivative of y(I) scaled by * 


10 

20 

30 

AO 

50 

60 

70 

30 

90 

100 

110 

120 

130 

lAO 

150 

160 

170 

130 

190 

200 

210 

220 

230 

240 

250 

260 

270 

280 

290 

300 

310 

320 

330 

340 

350 

360 

370 

380 

390 

400 

410 

420 

430 

440 

450 

460 

470 

480 

490 

500 

510 

520 

530 

540 

550 


94 



c* 

Co 

C« 

c* 

CO 

Co 

Co 

CO 

Co 

CO 

Co save 

Co CSAVt 

Co H 

CO 

Co 

CO 

Co 

CO 

Co 

c* 

Co HHIN 

c* 

Co 

C* 

Co HMAX 

C« EPS 

Co 

CO 

Co 

C* “F 

Co 

Co 

Co 

Co 

c* 

c* 

Co 

C* 

Co 

Co 

Co 

Co 

Co 

CO 

CO 

C* 

Co 

Co 

Co 

Co 

Co YMAX 

CO 

CO 

C* 

Co ERPOR 


Ho»j/faCTORIAU ( J) WHFRF H IS THE CURRENT 
STEP SIZE. ONLY Y(l.l) NEED BE PROVIDED BY 

the calling phorram on the first entry. 

IF IT IS desired To interpolate to non mesh points 
these values can be used. if the current step size 
is H and the value at t + e is needed* form 
s = E/H. ANO Then compute 

NO 

Y«I)(ToE) = SUM Y«J+l*n*Sooj 
J=0 

A BLOCK OF AT LEAST 13on FLOATING POINT LOCATIONS. 

N 03 floating point LOCATIONS USED BY THE SUBROUTINES. 
THE STEP SIZE TO d£ ATTEMPTED ON THE NEXT STEP. 

H MaY BE adjusted UP OR DOWN BY THE PROGRAM 
IN ORDER TO ACHEIVE AN ECONOMICAL INTEGRATION. 
HOWEVER. IF the H PROVIDED BY THE USER DOES 
NOT CAUSE A LARGER ERROR THAN REQUESTED* IT 
WILL BE USED. TO SAVE COMPUTER TIME* THE USER IS 
ADVISED TO USE A FAIRLY SMALL STEP FOR THE FIRST 
CALL. IT WILL be AUTOMATICALLY INCREASED LATER. 

THE minimum step size THAT WILL BE USED FOR THE 
INTEGRATION. NOTE THAT ON STARTING THIS MUST BE 
MUCH smaller than the average H EXPECTED SINCE 
A FIRST order method IS USED INITIALLY. 

THE maximum size TO WHICH THE STEP WILL 3E INCREASED 
THE ERROR TEST CONSTANT. SINGLE STEP ERROR ESTIMATES 
DIVIDED BY YMAX(I) MUST ^E LESS THAN THIS 
IN THE euclidean norm, THE STEP AND/OR ORDER IS 
ADJUSTED TO AChEIVE THIS. 

THE method INDICATOR. THE FOLLOWING ARE ALLOWED.. 

0 AN AOamS predictor CORRECTOR IS USED. 

1 A MUlTI-STEP METHOD SUITABLE FOR STIFF 

SYSTEMS IS USED. IT WILL ALSO WORK FOR 
NON STIFF SYSTEMS. HOWEVER THE USER 
must PROVIDE A SUBROUTINE PEDERV WHICH 

evaluates the partial derivatives of 

THE differential EQUATIONS WITH RESPECT 
TO THE Y»S. THIS IS DONE BY CALL 
pedervit.y.pw.m) . pw IS an n by n array 
WHICH Must be set to the partial of 

THE I-TH equation WITH RESPECT 
to the j dependent variable in PW(I*JJ. 
Pw IS actually stored in an H by m 
array where M is the value of N USED ON 
the first CALL TO this PROGRAM. 

E THE Same as case i* except that this 
subroutine COMPUTES THE PARTIAL 
derivatives by numerical differencing 
OF THE derivatives. HENCE PEDERV IS 
NOT CALLED. 

AN ARRAY OF N LOCATIONS WHICH CONTAINS THE MAXIMUM 
OF EACH Y SEEN SO FAR. IT SHOULD NORMALLY BE SET TO 
1 IN EACH COMPONENT BEFORE THE FIRST ENTRY, (SEE THE 

DESCRIPTION OF LPS.) 

AN array OF N elements WHICH CONTAINS THE ESTIMATED 
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580 
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C« 



ONE step error in each COMPONENT. 

« 

1110 

C« 

KFLAG 

4 

completion code with the following meanings.. 


1120 

C* 



+1 the step was succesful. 

« 

1130 

C« 



-1 THE step mas taken with H = HMIn, BUT THE 

« 

1140 

c« 



requested error was not achieved. 


1150 

C» 



-2 the Maximum order specified was found to 


1160 

c* 



BE TOO large. 


1170 

c* 



-3 corrector convergence COULU NOT BE 

« 

IIBO 

c* 



achieved for h .gt. hmin. 


1190 

c» 



_4 THE REQUESTED ERROR IS SMALLER THAN CAN 


1200 

C4 



BE handled for THIS PROBLEM. 

« 

1210 

c* 

JSTART 

4K 

input indicator with the FOLLOWING MEANINGS.. 

o 

1220 

c* 



-1 repeat the last Step with a new h 

« 

1230 

c* 



0 perform the first STEP. THE FIRST STEP 


1240 

c* 



must be done with this value of JSTART 

» 

1250 

c» 



SO THAT THE SUBROUTINE CAN INITIALIZE 


1260 

c« 



itself. 

« 

1270 

c« 



+1 take a new STEP CONTINUING FROM THE LAST. 


1230 

C4 



JSTART IS SET TO NO, THE CURRENT ORDER OF THE METHOD 

« 

1290 

c« 



AT FXIT. riQ IS ALSO THE ORDER OF THE MAXIMUM 


1300 

C4 



derivative available. 

« 

1310 

c* 

maxDER 

T^IE maximum derivative that should be used in The 


1320 

c* 



METHOD. SINCE THE ORDER IS EQUAL TO THE HIGHEST 

« 

1330 

c* 



DERIVATIVE USED. THIS RESTRICTS THE ORDER. IT MUST 


1340 

c* 



BE LESS than 13 FOR ADAMS AND 7 FOR STIFF METHODS. 


1350 

c* 

PM 

A 

BLOCK OF AT LEAST N«»2 FLOATING POINT LOCATIONS, 

•» 

1360 

c* 

IP 

A 

BLOCK OF AT LEAST N INTEGERS. 

« 

1370 


SUeROUTINE OirSflM (N«T.Y«SAVE:,CSAVe»H.HMlN.HMAX*£PS»MFf YMAXtEkROR, 1390 

IKFUAO.JSTARTtMAxOEPtPi^tlP) 1400 

OliMtNSION Y(13,l)« YmAX(1)» SAVE(13»D* ERROR(I)* PW ( 1 ) * A(13)» PE 1410 

1RTSTU2,r,3) , CSAYE(N*3)* IP(d 1420 

C* THE COthFiCTENTS IN PERTST ARE UEED IN SELECTING THE STEP AND * 1440 

C* ORDER* therefore ONLY ABOUT ONE PERCENT ACCURACY IS NEEDED. * 1450 

DATA PERTST/?.. 4.5. 7.333. 10. 42, 13.7, 17. 15. 1..1..1..1..1..U.2.. 12. 1470 

1. 24* » 37. p9, 53, 33, 70. 08.57. 97, 100. 9, 126. 7. 147. 3, 168. 8. 191. 4. 3. .6., 9 1480 

2.167.12.5.15.98.1.. 1..1..1..1..1..1..12..24..37.89.53.33.70.08.87. 1490 

397.100.9. 126.7. 147.3. 188.9.191.4.211. . 1.. 1 . . .5. . 1667. . 04 1 33 . . 0 0826 1500 

47,l.»l.,l.,l.,l.,l.,l.,l.,2.,l.,.3l57,.07407,.0139,.0021818,.00029 1510 

5. . 000035. .000 0037.. 00 00 0035/ 1520 

DATA A(2)/-1.o/ 1530 

IRET=1 1540 

KFLAG=1 1550 

IF lUSTflHT.LE.O) GO TO 50 1560 

C* BEGIN BY saving INFORMATION FOR POSSIBLE RESTARTS AND CHANGING * 1580 

c* H BY THC Factor R if the caller has changed H. all variables • 1590 

C* dependent on H must ALSO be changed. * 1600 

C* E IS A comparison for ERRORS OF THE CURRENT ORDER NO. EUP IS * 1610 

C» TO TEST FOP increasing T«E ORDER, EO<<N FOR DECREASING ThE ORDER. * 1620 

C* HNEW is the step size that has used on the last CALL. * 1630 

C »»«««« 4* *** 4 «•«>«»««»» «*«•<« »«»«««« 1640 

10 00 «:u 1 = 1. N 1650 
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DO 20 J=1.K 1660 

20 S>‘Ve(J,I)=Y(J,I) . 1670 

HOLD=HNfw 1690 

IF <tt.EQ.H0LD) DO TO <*Q 1690 

30 IREU=1 1700 

GO TO 750 1710 

40 NOOCUsNo 1720 

TOLD=T 1730 

IF <OSTftPT.GT.O) GO TO 300 1740 

GO TO 8n 1750 

So IF <'^STiPT,EQ.-l) GO TO 70 1760 

C aa «««« os 1770 

C» CM THE fiRsT call. THE ORDER IS SET TO 1 AND THE INITIAL * 1790 

C» DERIVATIVES ARE CALCULATED. * 1790 

Co»»«*o'»***-»»oh»*****-i»»o***o«««***«»*«o»******»'»**'»* 4'»*»'»*'»*4*****« ****** leOO 
RP=I.O 1910 

NQ=1 1920 

N3=''i 1930 

N4=N**2 1940 

call OIfFUN (N,T.Y.CSAVE(1.3) ) 1850 

00 60 1=1, N I860 

60 Y<«i.I)=CSAVE(I.3)*H 1870 

HNEW=H 1880 

K=2 1990 

GO TO lo 1900 

C *«»■»«****'*■'**■»« *«»*o **«•»«***■»<**«*««*«*•»**■*«*•»**•»«*•»**•»********■»* *■**•»■**** 1910 

C* repeat last step 0Y restoring saved INFORMATION. * 1920 

C*«*««<H»«****««(HHt-»*0*««««0'*44*»<H>*««*******'»«*'»**<HH»**0**0*****'»*«***»« 19 30 

70 IF ‘NQ.eO.NOOLD) JSTART=1 1P40 

T=TOLO 1950 

NQ=NU0LD 1960 

K=NO+l 1970 

GO TO 30 1980 

1990 

C* SET THt coefficients THAT DETERMINE THE ORDER AND THE METHOD * 2000 

c* TYPE, lheck For excessive order, the last two statements OF * aolo 

C» This StcTIoN SET IwEVAL .OT.O IF PW IS TO BE RE-EVALUATED * 20P0 

C» BECAUSE OF the order CHANGE, AND THEN REPEAT ThE INTEGRATION * 2030 

C» STEP IT IT has NOT yET B^EN DONE (IRET = 1) OR SKIP TO A FINAL * 2040 

C» scaling before FXIT if it HAS BEEN COMPLETED (IRET =2). * 2050 

C********”************#*** »***«■»*»■»*«*•»<*** *»***»*******■»**********•»**«*« 20&U 

80 IF (MF.EO.O) GO TO 90 2070 

IF InQ.pt.B) go to 100 2090 

GO 10 (230.240»250»260,270.280) . NO 2090 

90 IF iNO.fiT.l?) GO TO 100 2100 

GO TO (110, 120. 130. 140. 150, 160, 170. 180. 190, 200*210, 220) , NG 2110 

lOO KFLA6=-? 2120 

RETURN 2130 

««««««««««<»« «««««*««*««««««« «•«««««««« 2^4^ 

C» the following coefficients should be DEFINED TO THE MAXIMUM » 2150 

c» ACCURACY Permitted by the machine, they are. in the order used.. * 2i6o 

C» * 2170 

C* -1 * 2180 

C» -1/2. -1/2 * 2190 

C» -5/12, -3/4. -1/6 * 2200 
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C* -3/p«-U/1?,-1/'3,-i/ 24 * 22in 

C* -251/7<dU,-?S/24, -35/72, -5/48, -1/120 * 2220 

C* -95/280, -137/120, -<5/8, -ir/'^b, -1/40, -1/720 » 2230 

C* -19087/b04pQ, -49/40, -203/'i70, -49/102, -7/144, -7/1440, -1/5040 « 2240 

C* -5257/17280,-383/280,-469/540,-967/2980,-7/90,-23/2160,-1/1290, 2250 

C* -1/40320 2260 

C* -107001 '/ 36 28800,-76 1/560, -2953 1/30240, -267/640, -1069/9600,-3/ 160, 2270 

C* -13/6720,-1/8960,-1/362880 2280 

C* -4 165506/ 1451 520 0,-7 129/5040, -65 15/6040, -4523/90 72, -19/ 128, 2290 

C» —>013/1036800,-5/1344,-29/96768,-1/72576.-1/3628800 2300 

C« -134211865/479001600,-7381/5040.-177133/151200,-84095/145152, 2310 

C» -341693/1814400,-8591/207360.-7513/1209600,-121/193536, 2320 

C* -11/272160.-11/7257600,-1/39916800 2330 

C» -262747<;65/95800320 0,-83 71 1/55440, -190553/151 200, -341 747/518400, 2340 

C» -139381/604800,-2667907/47900160,-1903/201600,-10831/9676800, 2350 

C» -11/120960,-1/207369,-1/6652800,-1/479001600 2360 

C« * 2370 

C* » 2380 

C« -1 * 2390 

C* -2/3, -1/3 » 2400 

C« -12/25,-7/10,-1/5.-1/50 * 2410 

C* -120/274,-225/274,-85/274,-15/274,-1/274 * 2420 

C* -180/441, -<58/63, -15/36, -25/252, -3/252, -1/1764 4 2430 

110 A(l)=-1.0 2450 

GO TO 290 2460 

120 A(l)=-,50onnnooooo 2470 

a (3)=-0. 500000000 2480 

GO iO 290 2490 

130 A(l)=-0, 4166666666666687 2500 

4(3)=-0. 750000000 2510 

4(4)=-0. 1666666666666667 2520 

GO TO 290 2530 

140 a(l)=-0, 3750000000000 2540 

a (3) =-0 .9166666666666667 2550 

A(4)=-'1. 3333333333333333 2560 

a (5)=-0, 04166666666660687 2570 

GO TO 290 2580 

150 a<l)=-0. 3486111111111111 2590 

a<3)=-l. 0416666666666667 2600 

a<4)=-0. 4861111111111111 2610 

a (5)=-0. 1041666666666667 2620 

a(6)=-0. 008333333333333333 2630 

GO TO 290 2640 

IbO a(l)=-0. 3298611111111111 2650 

a(3)=-l, 1416666666666667 2660 

a(4J=-0. 625000000 2670 

a(5)=-0. 1770833333333333 2680 

a (6J =-0.0250000000 2690 

a (7)=-0. 001388888888888889 2700 

GO TO 290 2710 

170 a(l)=-0. 3155919312169312 2720 

a (3)=-l .225000000 2730 

a(4)=-0. 7518518519510519 2740 

a (5>=-0. 2552083333333333 2750 
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A(6)=-n.04af,lllllllllllll 2760 

A(7)=-0.o04«61111111111111 2770 

A(8)=-0.n001Q84126984l26<?84 2780 

GO TO 2qn 2790 

180 A (l)=-0. 9042245370370370 2800 

A (3)=-l. 296428571428485 2810 

A(4)=-0. 8685185185185185 2820 

A(5)=-0. 3357638888888888 2330 

A (6)=-0. 07777777777777777 2940 

A(7)=-0. 01064814814814815 2850 

A (8J=-0. 0007936507936507937 2860 

A {9)=-0. 0000248015873015873 2870 

GO TO 290 2380 

l90 A ( 1)=-O.2948f,8o0o4409l71 2890 

A (31=-1 .35892057142357 2900 

4 (4)=-0. 976554232804233 2910 

A (5)=-0. 41710750000 2920 

A (6)=-0. 1113541666666667 2930 

A (7)=-0. 0187500000 2940 

A (8)=-0. 001934523809523809 2950 

A(9)=-0. 000111607142857143 2960 

A ( 10 >=-0.000 002755731 92239859 2970 

GO TO 290 2980 

200 A ( 1 *=-0.2869754464285714 2990 

A<3>=-1, 41448412698413 3000 

A(4)=-l. 07721560846561 3010 

A (5>=-0, 498567019400353 3020 

4(6)=-0, 148437500000 3030 

A(7)=-n. 0290605709876543 3040 

A<8)=-0. 0037202380952381 3050 

A (9)=-0. 000299685846550847 3060 

A ( 10* =-0.0 00 01 37706596 11 9929 3070 

A (ip=-0.0 000 002755731 9223986 3080 

GO to 290 3090 

210 A< l>=-0. 280189564439367 3100 

A (3) =-l ,46448412698413 3110 

A(4)=-l. 17151455026455 3120 

A(5>=-0. 579358190035273 3130 

A (6>=-0. 188322861552028 3140 

A (7>=-0. 0414303626543210 3150 

A(8)=-0. 0062111447989418 3160 

A (9)=-0. 000625206679894180 3170 

A (10 >=-0.0000404 1740 1528512b 3130 

A < 1 1 >=-0.0000 0151565255731922 3190 

A (12> =-0.0 0 no 0 0025052 10838544 17 320 0 

GO lO 290 3210 

220 A(l>=-0. 274265540031599 3220 

A (3>=-l. 80993867243867 3230 

A(4>=-1. 26027116402116 3240 

A(5>=-0. 659234182098765 3250 

A (6>=-0.?3045800?645503 3260 

A(7)=-0. 0556972461052322 3270 

A (8>=-0. 00943948412698413 3280 

A(9>=-0. 00111927496693122 3290 

A ( 10> =-0.0000909391534391534 3300 
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A (11 )=-0.0 on (1048 2253 08641 9753 


3310 


A ( 12> =-0.0 00 non 150 3 126503 12650 


3320 


A ( 13) =-0.00000000200767569878681 


3330 


GO TO 290 


3340 

230 

A(l)=-1. 000000000 


3350 


GO TO 290 


3360 

2^0 

A { 1 )=-0. 6666666666666697 


3370 


A (3) =-0.1333333333333333 


3380 


GO TO 290 


3390 

250 

A <1 ) =-0.5454 545454545655 


3400 


A <3)=A (1 ) 


3410 


A (4) =-0.0909090909090 90 91 


3420 


GO TO 250 


3430 

260 

A (1 >=-0.480000000 


3440 


A (3)=-0. 700000000 


3450 


A (4)=-0. 200000000 


3460 


A (5)=-0. 020000000 


3470 


GO TO 290 


3480 

270 

A ( 1 >=-0.437956204379562 


3490 


A (3> =-0.821 16788321 16788 


3500 


A <4> =-0.3 1021 8978 1021898 


3510 


A (5) =-0.054 74452554744526 


3520 


A (6> =-0.0036496350364963504 


3530 


GO 10 290 


3540 

280 

A (I r=-0. 408 1632653061 225 


3550 


A (3>=-0. 9206149206349206 


3560 


A (4> =-0.4 166666666666667 


3570 


A (5> =-0.0992063492063492 


3580 


A (6) =-0.0 1190476 190476 19 


3590 


A (7> =-0.0 00566893424036282 


3600 

290 

K=NU+1 


3610 


IDOUO=K 


3620 


MTYH=(4-uF) /? 


3630 


ENQC=. 5/FLOAT (NO+1 ) 


3640 


ENQJ=.5/FL0AT (N(3+2) 


3650 


ENQA=0.5/FLnAT{NQ) 


3660 


PEPSM=EP5 


3670 


EUP=(PErtST (N0,MTYP»2> *PEPSH) *«2 


3680 


E=(PtRT5T (NO.MTYP. 1 > *PEPSH) ««2 


3690 


EDWN=(PirRTST (NQ.MTYP.3) OPEPSH) **2 


3700 


IF (tOWN.E(3.0) GO TO 780 


3710 


8ND=(EPs*EN03) **? 


3720 


IWEVal=mF 


3730 


GO TO (300.670) . IRET 


3740 


3750 

C» THIS Stl-TIOM COMPUTES TH£ PREDICTED VALUES BY EFFECTIVELY 


3760 

Co MULTIPUyINO the saved information by the pascal triangle 


3770 

Co MATRIX. 


3780 


3790 

300 

T=T+H 


3800 


00 310 J=2,K 


3810 


DU 310 J1=J,K 


3820 


J8=K-J1+J-1 


3830 


DU 3l0 1=1, N 


3840 

310 

Y(J2,I)=Y(J2,I)*Y(j2*1 .1 ) 


3850 


100 



C» UP TO 2 cnRPPCTOR iterations Ape TAKEM. convergence is tested DY « 3870 

Co REQUIRING TMe L2 NORM OF CHANGES TO 8t LESS THAN BNO WHICH IS » 3880 

C* OEPENUENt ON THE ERROR TEST CONSTANT. * 3890 

C« THt SlJM OF THE corrections I$ ACCUMULATED IN THE ARRAY * 3900 

c» FRROR(I). it is equal to the K-TH derivative of Y multiplied « 3910 

c* PY h»*k/ (FACTORIAL <K-1)*A (K) ) , aNO IS THEREFORE PROPORTIONAL * 3920 

Co TO the ACTUAL ERRORS TO THE LOWEST POWER OF H PRESENT. (H*oK) o 393O 

C O <0 -tto *0 O'* »***«****♦«« 394O 

00 Ji:0 T = 1 ,M ' 3950 

320 ERROR(r)=0.0 3960 

no AJO L=I'P 3970 

call niFFUM <N.T, Y.CSAVE ( 1 .3) ) 3980 

C*-»«**'»««**'»'»»*<t»****<nt»*oo**®«e****'»«**»****'»***'»*********'*»**'»'»**<‘-»*** 3990 

c* IF there has «E£n a change of order or there HAS been TROUBLE » 4000 

C* WITH CONvFRGENCE. PW is RE-EVAluAT£0 PRIOR TO STARTING THE * 4010 

C* corrector ITERATION IN THE CASE OF STIFF METHODS. I'*EVAL lb * 4020 

C* THEN SET TO -1 AS AN INDICATOR THAT IT HAS BEEN DONE. » 4030 

IE (IwEVAL.lt. 1) go TO 390 4050 

IT (MF.EQ.2) go to 360 4060 

call oEDERV (T.Y.PW*N3) 4070 

R=«(1)*H 4080 

DU 330 1=1. N4 4090 

330 HW<t)=PW(I)*R 4100 

C*******<>*^<M><M>'»*«* •>*•&*•»« **-«K>*^*0'»*<HH>'*'«H>‘«K>*****<HHHH>**** *■»♦*** *♦<>*<>**** 41^0 

c» ADO thc. identity matrix to the Jacobian and decompose into lu = pw * 4120 

C0««*«»*0«*4«0*4*»*»400»0«»<*«4«**«**.,.,*«0**»»*«(*«H»0*0-»*-»*******<»**«***«* 4 130 

340 DU 3bo 1=1. N 4140 

350 PW(I*(N3+l)-N3)=1.0't'Pw<I*(N3'Hl)-N3) 4150 

IWEVAl=-1 4160 

CALL NDIOIZ (N.N3.PW.IP) 4170 

IT (IP(N).NE.O) GO TO 390 4180 

GO TO 440 4190 

C» evaluate The JACORIAN into Pw by numerical differencing. R Ib THE * 4210 

C» CHANGE MADE TO THE ELEMENT OF Y. IT IS EPS RELATIVE TO Y WITH * 4220 

C* A MINIMUM OF EP5**2. F STORES THE UNCHANGED VALUE OF Y. * 4230 

C*04-»*«4***»****«*»OC*«***o4*A»*o4«O **«***»*■»*•»»«*»*■»**■»■»***■»■»** ******** 42 AO 

360 DO 380 J=1.N A250 

T=Y(1.J) 4260 

H=Eps»AMAX1 (EPS.ABS (F) ) 4270 

Y ( 1 . J) =Y (1 . J) +R A280 

D=A ( 1 ) oh/r 4290 

LALL DIFfuN (N.T.Y.CSAVE(I.D) A300 

DO 370 1=1. N A310 

370 PW(I•^(J-1)»N3)=(CSAVE(I.1)-CSAVE(I.3) )»0 4320 

380 Y(1,j)=f 4330 

GO TO 340 4340 

390 DO 40n 1=1. N 4350 

400 CSAve(i,i)=y(2.I)-CSAVEU.3)*H 4360 

IT (MF.EO.O) GO TO AlO 4370 

CALL M^IOPr (N.N3 .PW.LSAvE<1,1) .IP) 4380 

C* CORRECT and compare DEL. THE L2 NORM OF ChANGE/YMAX. WITH 0NO. * 4400 
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C» ESTIMATt ThF value OF THE L2 NORM OF THE NEXT CORRECTION BY * 4A10 

C« BR»3*DtL AkjO COMPARE WITH RND. IF EITHER IS LESS* THE CORRECTOR * 4420 

C» IS SAI^ TO have converged. « 4430 

4l0 DtL=0,0 4450 

OU 42n 1=1, N 4460 

T ( 1 ,T)=Y (I , I) *A ( I) *CSAVE (I ,1 ) 4470 

T (2, I )=Y (2, I ) -CSAVE ( I , 1 ) 4480 

tRROo( I)=ERROR ( I ) ■•■CSAVEI 1, 1 ) 440 0 

DEL=DEL+(CS4VE(I»l)^TMAX(n)**2 4500 

420 LONTINUP • 4510 

It (U.RE.2) BR=AMAX1 « .0*BR,OEL/OEL1 ) 4520 

OtLl=nEL 4530 

IF (AminI (OEL,BR*OEt-*2.0) .LE.BNO) GO TO 480 4540 

430 CONTlMirE 4550 

C» the corrector ITERATION FAILED TO CONVERGE IN 2 TRIES. VARIOUS * 4570 

C* possibilities are checked for. if H IS ALREADY HMIN AND * 4580 

C» THIS IS either ADAMS METHOD OR THE STIFF METHOD IN WHICH THE * 4590 

c» MATRIX RW has ALREADY BEEN RE-EVALUATED, A NO CONVERGENCE EXIT * 4600 

C» IS taken, otherwise the matrix PW is re-evaluated and/or the * 4610 

Co step is reduced to try and get Convergence. * 462o 

440 T=TOLO 4F.40 

IF (aHS(H) .LE.HMIN*!. 00001. AND. ( IWEVAL-MTYP) .LT.-l) GO TO 450 4650 

IF I tMF.FO.O) .OR. (IWEVA l.NE.O) ) H=H*0.25 4660 

IwEVAL=mf 4670 

IRET1=2 4680 

GO TO 7sn 4690 

450 KFLAG=-3 4700 

NQ=NwOLD 4710 

460 DO AIO 1=1, M 4720 

DU 470 J=1,K 4730 

470 YlJ.I)=SAVE(J,I) 4740 

h=HOLD 4750 

JSTaHT=NO 4760 

return 4770 

«»«»»««««»««■»»««»«««« 4780 

C» THE CORRECTOR CONVERGED AnD CONTROL IS PASSED TO STATEMENT 520 » 4790 

C* IF THE error TFST is O.K., AND TO 540 OTHERWISE. * 4800 

C« IF the step IS O.K. IT IS accepted. IF IDOUB HAS BEEN REDUCED » 4310 

C* TO ONE* A test is made to see if the STEP CAN BE INCREASED * 4820 

C» AT THE CURRENT ORDER OR BY GOING TO ONE HIGHER OR ONE LOWER. * 4830 

c« SUCH A Change is only made if the step can be increaseo by at * 4840 

C* least i.l. IF MO CHANGE IS POSSIBLE IDOUB IS SET TO 8 TO • 4850 

C* PREVENT FUtHER TESTING FOR 8 STEPS. * 4860 

C» IF A change is POSSIBLE. IT IS MADE AND lOOUR IS SET TO * 4870 

Co NO + 1 TO PPFVENT further tESING FOR THAT NUMBER OF STEPS. o 4930 

C* IF THE ERROR WAS TOO LARGE, THE OPTIMUM STEP SIZE FOR THIS OR * 4890 

Co lower order is COMPUTED, AND THE STEP RETRIED. IF IT SHOULD o 49OO 

Co fail twice MORE IT IS AN INDICATION THAT THE DERIVATIVES THAT o 4910 

Co have ACLUMulaTFO in THE Y ARRAY HAVE ERRORS OF THE WRONG ORDER * 4Q?0 

C* SO THE FIRST DERIVATIVES ARE RECOMPUTED AND THE ORDER IS SET o 4930 

Co TO 1. o 4940 

Co ««»«««««»««««»«««« «««««««««<» 4950 
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460 0=0.0 4Q60 

00 •*'^0 t = l *.v 4Q70 

490 0=0+ (FOOOO ( I) /Y’^ftX ( I ) > »»2 4C)80 

I9£VAl. = n 4990 

IF (O.Gt.E) r,o TO 530 5000 

IF (O.LT.3) 00 TO 510 5010 

c» COvPteTc Tk-e CoPRFCTlON OF the higher O^DER derivatives after a * 5030 

C* SOCCEStUL STEP. » 5040 

DO RUO j=3,K 5060 

DO 500 1=1 5070 
500 Y IJ, I ) =Y ( J. I ) +A ( J) “tRROR ( I ) 5060 

51 I) KFLAG = *1 5090 

H.yEW=h 5100 

IF UDouP.LE.n GO TO 5110 

IOOub=Ipou6-1 5130 

IF (lOOllR.GT.l) GO TO 990 5130 

DO =90 t=l«M 5140 

530 C5AVE(T.2)=ERR0R(I) 5150 

GO TO son 5160 

C» REDUCE Tm£ FAILURE FLAG COlj^T TO CHECK FOR MULTIPLE FAILURES. * 5180 

C+ RESTORt T TO ITS ORIGINAL VALUE AND TRY AGAIN UNLESS THERE HAVE » 5190 

C® three t-AlLiiPES. In THAT CASE THF DERIVATIVES ARE ASSUMED TO HAVE » 5200 

C® accumulated errors so A RESTART FROM THE CURRENT VALUES OF Y 15 ® 5210 

C® tried. THIg is CONTINUED UNTIL SUCCESS OR H = HMIN, * 5220 

530 KFLAG=KFLAG-2 5240 

IF lAHS(H) ,LF. (HmIn*1.00OQ1) ) GO TO 740 5250 

T=TOLO 5260 

IF (KFLAG.LF .-5) GO TO 710 5?70 

C® PHlt PR2. and Pr3 will CONTAIN THE AMOUNTS 5Y WHICH THE STEP SIZE ® 5290 

C* CHOULO BE DIVIDFO AT ORDER OnE LOWER* AT THIS ORDER. AND AT ORDER « 5300 

C® ONE HIUHFp rESPFCTIVELY. * 5310 

540 PR2= (D/F) ®*ENQ2*i .2 5330 

PR3=1.E*20 5340 

IF < (NQ.o-E.MAXDER) .OR. (KFLAG.LE.-1> > GO TO 560 5350 

D=0.0 5360 

DO =50 1=1 ,M 5370 

550 0=0+ ( (FRROR ( 1 1 -CSAVE (I *2) ) /ymAX ( I ) » ®*2 5380 

PR3= (O/FijP) ®®EN03®1 .4 5390 

560 PR1=1.£*20 5400 

IF (NO.lp.I) GO TO 500 5410 

D=0.U 5420 

DO 9^0 1=1, M 5430 

570 0=U* ( V (K, I ) /YMAX ( I ) )®«2 5440 

PR1= (O/FpwN) »®E n01 »1 .3 5+50 

560 CONliNUF 5460 

IF IPR2.LE.PD3) GO TO 640 5470 

IF (PP3.LT.PR1) GO TO 650 5480 

590 P=1 .O/Aujxi (BPl , 1 .E-4) p490 

NEWU=NQ-I 5500 
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600 lOOUb=H 5510 

IF ( iKFLAG.Frv. 1 ) .and. (R.LT. ( 1. 1) ) ) GO TO 6<?0 5550 

IF lNbWO.LE.NO) GO TO 620 5530 

C* COMPuTt ONF AOniTIONAL scaled DtRIVATIVE IF ORDER IS INCREASED. * 5550 

DO biO 1 = 1. N 5570 

610 r<^E«0*l.I)=FPROR(I)*A(K)/FLOAT(K) 5550 

6i:0 K=NEW3*i 5590 

IF tKFLAG.EO.l) GO TO 660 5600 

H=H*R 5510 

IREU = 3 5620 

GO TO 790 5630 

630 IF lNEwn.EO.MO) GO TO 300 5640 

NQsNtwQ 5650 

GO TO Go 5660 

640 IF ‘PR2.RT.PP1) GO TO 590 5670 

NEWO=NQ 5680 

r=1.0/AM4X1 (PR?,1.E_4) 5690 

GO TO 300 5700 

650 R=1.U/AM4X1 {OP3.1.E-4) 5710 

NEWO=NQ+i 5720 

GO TO 600 5730 

660 I»Ei=2 5740 

R=AMIN1 (R.HMAX/ABS (H) ) 5750 

H=H*R 5760 

H.NE"=R 5770 

IF INQ.En.NFwQ) GO TO 670 5780 

NQ=NEwO 5790 

GO TO 80 5800 

670 Rl=i.O 5010 

DO 660 J=2,K 5920 

Rl=Rl»9 5930 

DO 680 1=1, N 5840 

680 Y lo. I )=Y ( J, I) *Rl 5850 

I00U6=K 5060 

690 DO /UO 1=1, N 5870 

700 YMAX (I)=4M4X1 ( YMAX (I ) .A0S (Y(l ,1) ) ) 5880 

jSTART=N0 5890 

RETURN 5900 

710 IF (NQ.rt.D (50 TO 720 5910 

IF ‘ABS (H) .LE.2.*hMIN) GO TO 780 5920 

GO ro 540 5930 

720 R=H/hOLD 5940 

DO 730 1=1, N 5950 

Ytl,I>=SAVF(l«I) 5960 

730 Y‘<i,I)=SAvE(2,I)*9 5970 

NO=l 5980 

KFLAU=1 5990 

60 TO 80 6000 

740 KFLAUs-1 6010 

HN£W=H 6020 

JST8RT=N0 6030 

return 6040 
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c« THIS StCTION SCALES ALL VARIABLES CONNECTED WITH H AND RETURNS * 6060 

C« TC THE cNTrsiNo SECTION. * 6070 

760 H=S1'^N(AmAX1 (Hmin*AMIN1 'AQs(m) .hmAX) ) tH) 6090 

Rl=i*0 6100 

00 T60 j=2.k 6110 

R=H/Hf)LO 6120 

Rl=Rl*o 6130 

00 760 1=1 .N 6140 

760 Y(J.I)=SAVF(J,I)*R1 6150 

no ItO 1 = 1. Kj 6160 

770 Ytl,I)=SAVF(l.l) 6170 

I00UO=K 6180 

r,0 TO (40.300.630). IHETI 6190 

780 KFLA0=-4 6200 

GO TO 460 6210 

END 6220 

SUBHOUTinE MPIOlZ (N.NDIM.A. IP) 6230 

— HOLER'S "OFCOMP". 6240 

6250 

matrix Tpi ANGULAPIZATION by GAuSSIAN elimination. 6260 

6270 

INPUT... 6280 

N = OpOro OF MATRIX 6290 

NOIM = DECLARED DIMENSION OF ARRAY A, 6300 

A = matrix to re TRIANGuLARIZED. (FOR STIFF METHODS, A IS SINGLE 6310 

PRECISION) ALL other VARIABLES APE DOUBLE PRECISION.) 6320 

OUTPUT... 6330 

A(I.J). I.LF.J = UPPER triangular FACTOR, U. 6340 

A(I,0), I.GT.J = multipliers = LOWER TRIANGULAR FACTOR, I-L. 6350 

IP(6). K.LT.N = INDEX OF K-TH pivot ROW. 6360 

IP(N» s (-1)4*{NUMBEH of INTERCHANGES) OR 0. 6370 

USE 'SOLVf, (NDI02Z) TO OBTAIN SOLUTION OF LINEAR SYSTEM. 6380 

OETERM(A) = IP (N) »A ( 1 , 1) *A (2,2) • , . .»A (N.N) . 6390 

IF IP(N) = 0. A IS SINGULAR, SOLVE WILL DIVIDE BY ZERO, 6400 

6410 

DIMENSION A (NDIM.NOIM) , IP(NDIM) 6420 

IP(N)=l 6430 

DO ou K=|,N 6440 

IF (K.fO.m) go to 50 6450 

KHl=K*l 6460 

M=G 6470 

DU 10 I=KP1,N 6480 

IF ( ABS ( A ( I ,K) ) .GT. ABS ( A (H.K) ) ) M=I 6490 

10 UONTINUE 6500 

IH(K)=m 6510 

It- (M.nE.K) IP(N)=-IP(N) 6520 

T=«(M.k) 6530 

A IM.rC) =A (k.K) 6540 

A(K,K)=T 6550 

It- (T.pq.O) go to 50 6560 

DO 20 I=KP1.N 6570 

20 A(I,k)=-A(I,K)/T 6580 

DO 40 J=KP1.N 6590 

I=A (M,U) 6600 
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“ =A (K* J) 6610 

**(K.j)3T 6620 

IP" (T.EO.O.) GO TO 40 66T0 

00 30 I=KP1,N 6640 

30 4(t.j)-A(I,J)+A(I,K)»T 6650 

40 OONtINUP 6660 

50 IP (4 (K,K) .EQ.O. ) lP(N)s0 6670 

60 CUNTImuE 6680 

RETUKN 6690 

ENO 6700 

SUBMOUTt^E NOIOPZ (N,N0IM,a,8,ip) 6710 

— f<CLEP'S "SOLVE". 6720 

6730 

solution of linear system* 4»X = 8. 6740 

6750 

Input*.. 6760 

N = OMOFO of matrix. 6770 

NOIM = OECLAPEO nIMENSION OF ARRAY A. 6780 

A = TrIanGULARIZEO matrix obtained from »DEC0MP« (NDIOIZ). 6790 

P = right Hand side vECTOr. 6800 

IP = PIVOT VfCTOR obtained from iOECoMP*. 6810 

OUTPUT... 6820 

8 = SOLUTION VECTOR. X. 6830 

DIMENSION A (NDIm.NDIM) » 8(NDlM). IP(NOIM) 6840 

IF (N.EO.l) GO TO 30 6850 

NM1=N-1 6860 

00 iU K=i,mm1 6870 

KPI=K+1 6880 

M=IP(K) 6890 

T=B(M) 6900 

B(M)=R(K) 6910 

H(K)=t 6920 

DU 10 i=KPl.N 6930 

10 8 ( I )=a ( I) +A ( I ,K) »T 6940 

DO 80 KR=1,nm1 6950 

KM1=ai_kB 6960 

K=KM1*i 6970 

B lN)=q (K) /A (K.K) 6980 

T=-8(K) 6990 

DU 20 1=1, KMl 7000 

20 Ba)=R(I)*A(I,K)*T 7010 

30 B(1)=R(1)/A(1,1) 7020 

RETUkN 7030 

END 7040 
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